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Abstract 


A general technique for solving Maxwell’s equations exactly, based on expansion 
of the solution in a complete set of vector basis functions has been developed. These 
vector eigenfunctions are derived from the complete set of separable solutions to the 
scalar Helmholtz equation in a particular coordinate system and are shown to form a 
complete set. The method is applicable to a variety of problems including the study 
of near and far field electromagnetic scattering from particles with arbitrary shapes, 
plasmon resonances in spherical nanoparticles with spherically concentric ‘shells’ and 
the calculation of plasmon resonances in the sphere-plane geometry. An exact method 
for solving the inhomogenous Maxwell’s equations (i.e., in the presense of charges and 


currents) is also outlined. 
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Chapter 1 


Introduction 


The subject of obtaining the solution of Maxwell’s equations in a boundary value 
problem has kept scientists occupied for more than one hundred years [25]. It is one 
of the few ‘classic’ problems that continue to attract the attention of researchers. 
Many contemporary problems of interest involve our knowledge of the exact solution 
to these equations. However, the complexity of these equations has defied analytic 
solutions for all but the simplest cases. 

To state the problem more correctly, one needs to understand the length-scales 
involved. In a scattering problem for example there are two lengths involved: d, the 
physical size of the scattering object and \, the wavelength of electromagnetic waves 
in question. When d > A, the problem is readily solved. The method of solution, 
called geometrical optics, has been known for at least a few hundred years. In fact 
we do not even need to invoke Maxwell’s equations directly. The other limiting case 
is when d < A. This regime is called the quasistatic limit or nano-optics. Although 
mathematically involved, the solution in the quasistatic limit is obtained by solving 
Laplace’s equation. The problem is far more complicated when d ~ X. It is in this 
regime that one is required to solve Maxwell’s equations exactly. 

The problem involves the solution of a set of partial differential equations (The 
Maxwell Equations) for the coupled vector fields for the electric field E and the mag- 
netic field H. Although the general theoretical aspects of these equations are fairly 
well understood at the moment [25] [30], there are still some formidable mathemat- 
ical or computational problems to overcome in the area of solving the time-varying 
boundary value problem. Whereas the problem of acoustic (scalar fields) diffraction 
has been solved in a number of cases such as the strip, elliptic cylinder, hyperbolic 
cylinder, wedge, prolate and oblate spheroids, etc. [13], the corresponding solution 
for the vector fields are still by and large unconquered. Often the solution is worked 
out only in the quasistatic approximation by solving Laplace’s equation (homogenous 


case) [4] [3]. 


1.1 Historical Development 


In the absence of sources, these coupled partial differential equations have been solved 
exactly, analytically, only for a few isolated cases. The method of solution involved 
in these cases requires the knowledge of certain ingeneous transformations, suitable 
only for the problem in question. Usually these are dictated by the simple geometric 
shapes of the boundaries in question. Often these boundaries have to correspond to 
the coordinate surfaces in order for the boundary conditions to be separable. The 
simplest canonical case of scattering of plane waves reflecting/refracting at a planar 
interface [22] is commonly known by the name of Fresnel’s Equations. The most 
celebrated analytic solution to the scattering problem is that of Mie scattering, the 
scattering of electromagnetic waves from a dielectric sphere [12], originally due to G. 
Mie (1908). 

Immediately after Mie obtained the general solution to the sphere problem, Debye 
was able to formulate the same problem in terms of a pair of coupled scalar func- 
tions [27], in which the electric and magnetic fields were expressed. The method 
was again suitable for boundary value problems with spherical boundaries only. 
Hanzen [19] [20] [18] developed a technique for addressing the problem of radiation 
from antennas using a special type of transformation. The subject was developed fur- 
ther by Stratton [38] who solved the problem of Mie scattering using a set of vector 
functions derived from the solutions of the scalar Helmholtz equation. Subsequently, 
Aden and Kerker [1] used the formalism to solve the problem of electromagnetic scat- 
tering from two concentric spheres. Although Stratton [38] had made considerable 
contributions to this area of basis function expansion, the method of solution was still 
considered as just another ‘more elegant’ way to solve the Mie problem, analogous 
to the approach of Debye. The crucial idea that was missing at this stage was the 
notion of mathematical completeness of basis function expansions. Furthermore, the 
algebraic difficulties one had to overcome restricted the Stratton approach to solv- 
ing the spherical scatterer problem with incident light approaching along the z-axis 
only. Progress in this subject was negligible beyond this stage and the problem is 
reproduced in its original form even in contemporary text-books [40] [11]. 

The problem of diffraction of a plane wave at normal incidence on a circular cylin- 
der was originally solved by Rayleigh and since then his solution has been generalized 
and extended to plane waves at oblique incidence [41]. The method of basis function 


expansion was again used in this problem. Similarly, the problems of electromagnetic 


and acoustic scattering from a semi-infinite cone [36] and a semi-infinite body of rev- 
olution [35] were obtained by the method of basis function expansion. For the special 
case of the paraboloid, the solutions were ‘exact’. The generality of these approaches 
were not apparently proved (or even appreciated) beyond the specific geometries these 
earlier researchers were interested in. 

In the mid 1960’s there began a growing interest in the study of Maxwell’s equa- 
tions as a purely mathematical problem [25] [30]. Perhaps the growing use of so- 
phisticated electromagnetic devices testified beyond doubt the correctness of these 
equations for macroscopic electromagnetic phenomena. The theory of Maxwell’s 
equations was now to go through the rigors of verifying its mathematical founda- 
tions. The mathematical apparatus of topology and modern analysis was available 
to study the mathematically interesting properties of these equations. For example, 
it was found that any electromagnetic response within a perfectly conducting cavity 
(resonator) could be expressed in a complete set of functions (that is, the existence 
of such functions was proved) [25] [30]. 

The Sphere-Plane model consisting of a dielectric sphere close to a semi-infinite 
half-space substrate has been used to model the effect of irregularities in a surface 
that can lead to high local field enhancements [3] [4] [33] [34] [39]. The problem has 
been solved in the quasistatic approximation except in the work of Takemori [39] who 
obtained the exact solution using a Greens function approach with higher orders of 
scattering from the sphere and the plane. 

More recently, primarily due to the availability of increasingly powerful comput- 
ers, our approach for solving the Maxwell equations have taken a rapid ‘numerical’ 
path. It is all too common to try to solve a boundary value problem entirely numeri- 
cally. One uses finite difference equation approximations of the Maxwell equations on 
a (often non-linear) grid within a defined region of interest. Although such techniques 
are extremely powerful for solving problems with absolutely no symmetry, they are 
often very resource hungry in terms of memory(core) and computation time. Often 
researchers have to resort to solving the two-dimensional analogs (or some approxi- 
mation of it) to try to understand the physics in a specific situation. 

In the last decade there has been a rapid development in the area of scanning probe 
microscopies. As experiments have become more sophisticated, including the coupling 
of light at the active junction in a scanning probe microscope, considerable efforts 
have been directed towards modelling an illuminated tip-substrate geometry [14] [28] 
[10] [9] [8] [23]. Whereas the solution in the quasistatic limit is easy to solve [14], 


it is difficult to justify. The use of complex dielectric functions and the fact that 
the tip and substrate boundaries (the scattering objects in this problem) extend to 
infinity (so the condition that the dimensions are much smaller than the wavelength 
of the optical field in question does not hold true) are not totally justified. Attempts 
are currently underway to obtain ‘exact’ solutions to these problems. At the present 
moment, the available techniques to solve this problem too have remained by and 
large unsatisfactory. 

Some other notable geometries for which the exact solution of electromagnetic 
scattering have been obtained in three dimensions include the prolate and oblate el- 
lipsoids. Certain two dimensional cases, such as the circular, parabolic and elliptic 
cylinders, are also amenable for exact solutions [13]. However, for arbitrary geomet- 
rical shapes one has to resort to solving the Maxwell equations by finite element 
numerical techniques. 

Finally, the problem of solving the inhomogenous Maxwell equations in general 
have seen very little progress. The usual approach has been to reformulate the prob- 
lem in terms of the scalar and the vector potentials in a given gauge. One uses Greens 
function techniques to solve for the vector potential, given the current densities. The 
scalar potential is similarly solved from the knowledge of the charge densities. The 
mathematical apparatus is not sufficiently developed in this area and the problem is 


by and large unsolved, exactly. 


1.2 Outline of The Thesis 


In Chapter 2 we introduce the notion of the vector basis functions and prove that 
they form a complete set in a given coordinate system. The proof on completeness is 
followed by a discussion on the relevant boundary conditions in a scattering problem. 
In Chapter 3 we go into the detailed algebra of trying to expand an arbitrary plane 
wave in the basis set. Only the key results are highlighted and the bulk of the 
algebra is carried out in the appendices. The chapter on the Shell-Problem is the first 
application we consider in which the basis functions are used to obtain the solution. 
Comparison of the calculations are made with several experiments that are published 
in the literature. The next chapter on ‘Scattering from non-spherical objects’ uses the 
same basis function technique to solve the problem of scattering from an elongated 
object whose boundaries do not conform to the coordinate surfaces. Previously, this 


kind of calculation could be done by using completely numerical techniques only. We 


also look into the solution of the sphere-plane model and compare our calculations 
with experiments reported in the literature. Finally we conclude by highlighting the 
future directions in which this research could continue to establish itself as a powerful 


general technique for solving Maxwell’s equations. 


Chapter 2 


General Solution of Maxwell’s Equations 


2.1 Introduction 


Many macroscopic electromagnetic phenomena have been successfully described by 
Maxwell’s Equations. Over the past one hundred years, it has been put to test 
through the rigors of scientific scrutiny. Strictly speaking, the theory is applicable in 
a macroscopic sense, where one can assume the knowledge of ‘bulk’ electromagnetic 
properties of materials such as the dielectric function. However, the theory has been 
so successful that scientist have often taken the liberty of extrapolating the region of 
applicability to the nanometer and sub-nanometer scales. 

The electromagnetic field at a time t, and at any point r in a medium character- 
ized by a dielectric function ¢(r,t) and magnetic permeability p(r,t) is described by 
Maxwell’s Equations. The differential form of the equations in terms of the electric 
field E, the magnetic field H, the electric displacement D, the magnetic induction B, 
and the sources: J, the current density and p, the charge density, are defined below, 


for notational consistency within this thesis: 


V -D(r,t) = (r,t) (Gauss’s Law) (2.1) 

V x E(r,t) = ee (Faraday’s Law) (2.2) 

V -Bi(r,t) =0 (No Magnetic Monopoles) (2.3) 
OD(r, t) 


V x H(r,t) = J(r,t) + (Ampere’s Law) (2.4) 


Ot 

These four equations, along with the constitutive relations: B = wH and D = cE 
are in principle capable of describing all macroscopic electromagnetic phenomena. 
In practice, given an electromagnetic phenomenon describable as a boundary value 


problem, it is a formidable task to solve these equations exactly. At present, only 


finite element numerical techniques are available to solve a general problem. These 
numerical calculations are sufficiently resource-hungry in terms of memory and total 
computation time that in most cases only the simplest geometries and often only the 


two-dimensional analogs are attempted for a solution. 


2.2 Linear Homogenous Isotropic Medium 


For the interesting subclass of problems that are related to the scattering of electro- 
magnetic radiation by dielectric and/or magnetic boundaries, the general equations 
assume a much more symmetrical form. Within this approximation of a linear, time- 
invariant, homogenous and isotropic medium without any external sources, Maxwell’s 


equations assume the form: 


V-E(r,t) = 0 (2.5) 

V x E(r,t) =— a (2.6) 

V -B(r,t) =0 (2.7) 
AE(r, t) 


V x H(r, t) = cE(r,t) +€ (2.8) 


Ot 
Where we have assumed the possibility of an induced current density in the medium 
equal to oE, where o is the conductivity of the medium under consideration. Also, 
for harmonic time-dependence of the external fields, e~“*, where w is the angular 
frequency, the time derivatives could be factored out, and we obtain the set of cou- 


pled partial differential equations in the spatial coordinates alone, with the angular 


frequency of the excitation as a parameter : 


V-E(r,w) =0 (2.9) 
V x E(r,w) = twuH(r, w) (2.10) 
V - Bir,w) =0 (2.11) 


V x H(r,w) = oE(r, w) — weE(r, w) (2:12) 


From now on, we shall use a simplified notation by suppressing the argument list 


(r,w). Thus Equation 2.12 becomes: 


V x H+ iw(e+ )E=0 (2.13) 
w 
And from Equation 2.10 we obtain: 
1 
H = —VxE (2.14) 
Wf 


Substituting the expression for H from Equation 2.14 into Equation 2.13, and in- 
troducing the complex dielectric function é¢ = ¢« + %, we obtain the Vector 


Helmholtz Equation for the electric field vector: 
VxV x E-w*péE=0 (2.15) 


Similarly, taking the curl of both sides of Equation 2.13 and substituting for V x E 
from Equation 2.14, we obtain the corresponding Vector Helmholtz Equation for the 


magnetic field intensity vector: 
VxVxH-w*péH = 0 (2.16) 


These equations are totally symmetrical in the field variables and they also decouple 
the electric and magnetic fields. It is the solution to the Vector Helmholtz Equations 
for specified boundary and radiation conditions that describes the scattering of elec- 
tromagnetic waves. However, these equations are vector partial differential equations 
which are sufficiently difficult to solve in general. The most well known non-trivial 
problem that has been solved exactly, analytically, is that of scattering from a sphere 
(Mie Scattering, 1908). It is the purpose of this thesis to develop the general tech- 
nique of vector basis function expansion, for the solution of both the near and far 


fields in electromagnetic scattering from finite sized objects. 


2.3. Solution Space of Vector Helmholtz Equation 


The theory of solving the Scalar Helmholtz Equation, also known as the Wave 
Equation, is a very well developed subject [29] [5] [43] . This equation describes the 


propagation of scalar waves, such as acoustic waves, in a medium. 


Vw(r,t) + v(r,t) = 0 (2.17) 


When expressed in certain orthogonal curvilinear coordinate systems, such as the 
cartesian coordinate system or the spherical polar coordinate system, under assump- 
tion of sinusoidal time-dependence, the solution 2)(r,t) can be obtained by the tech- 
nique of separation of variables. The separation procedure reduces the partial dif- 
ferential equation to several ordinary differential equations. The separated equations 
can often be cast in the form of the well known Sturm-Liouville eigenvalue problem for 
second order ordinary differential equations, so that the solution space is guaranteed 
to be complete for the scalar Helmholtz equation. 

For the moment let us assume that we have obtained a complete set of scalar 
functions {7} that are solutions to the scalar wave equation (Eqn 2.17). Let us 


introduce the diffraction equation : 
V(V:-G)-VxVxG+hk?G =0 (2.18) 


where k? = w?é. The diffraction equation is satisfied by the electric field E and the 
magnetic field H that satisfies the Vector Helmholtz equations since the ‘extra’ term 
V -G will be zero for fields that have zero divergence. So, Equation 2.18 is consistent 
with Maxwell’s equations for electromagnetic waves. Let us define a vector function 


that is obtained by taking the gradient of the scalar function ~: 

L=V~w (2.19) 
L satisfies the diffraction equation if y is a solution to Equation 2.17: 
V(V-L)-VxVxL+hk?L 
V(V-Vb)-Vx Vx Vb +t kh?Ve 
V(W7b +k?) ( Since V x Veb = 0) 
0 (Since V2a + k?y = 0) (2.20) 

(221) 


Now let us assume that there exists a vector function M, with zero divergence, 
V:-M=0, that is a solution to the diffraction equation(Eqn 2.18). Let us con- 
sider the vector N = 7V x M. Assuming k is a constant(homogenous medium), we 


obtain: 


V(V-N)-VxVxN+k?N 


1 
= giV(V VxM)-VxVxVxM + kV x M) 
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= Vx(-V x(V x M)+k’M) (Since V-(V x M) = 0) 
= Vx(V(V-M)—Vx(VxM)+k?M) (Since V-M =0) 
= 0 (Since M satisfies Eqn 2.18) 222) 


Thus we show that N, as defined above, also satisfies the diffraction equation. Also 
V-N = 0, since divergence of curl is zero. So, if we had started with postulating 
the existence of a vector function N with zero divergence, we could show that there 


exists M with zero divergence. Thus from symmetry arguments alone, we could write 


Ni = cV x N. Specifically, since N = iV x M, since & is a constant (by assumption 


of homogeneity of the medium), 


1 

VN: = EYXVxM (2.23) 
1 

Von = peM (Follows from Eqn 2.18) (2.24) 


Thus M = Vv x N. Clearly, M and N are distinct from L, since the latter has non- 
zero divergence in general. So L must be linearly independent of {M,N}. It is up to 
us to create a vector function M(or equivalently N) from the given scalar function ¥, 
such that it will have zero divergence and will satisfy the diffraction equation. The 
point to note is that we could arrive at more than one set of functions {M,N}, and 
it is the relative algebraic convenience that will dictate the choice of a particular set. 

As a concrete example, let us consider M = V x ay, where y is the given scalar 
solution and a is an arbitrary constant vector. The divergence condition is satisfied, 
since divergence of curl is identically zero: V-M = V-(V x ay) = 0. Using 
the operator identity (the author has verified this ‘identity’ for the spherical polar 
coordinate sytem and the cylindrical coordinate system in addition to the rectangular 


coordinate system), V(V-M)—V x V x M= V?M, we can write: 
V(V-M)-VxVxM+k*?M 
= V’M+ih’M 
= V7(V x av)+ k?(V x av) 
= V(Vdxatv(V xa)) +k? (Ved xatv(V xa)) 
= V*(Vw xa)t+kh?(Vu xa) (Since V x a= 0) 
= (V7(Vedb)+h?Vod) xa (Since V? does not act on a) 
= (V(V7p +R) xa 
=f i) (Since x satisfies Eqn 2.17) (2:25) 
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Also, M=V xawv=Vwxa=Lxa. So,M-L=0,1.e. L and M are orthogonal. 
Thus, given a countably infinite set of particular solutions to Eqn 2.17, {v,}, that 
are finite, continuous, single-valued and with continuous partial derivatives; associ- 
ated with each vy, one can obtain a triplet of mutually non-coplanar vector solutions 
{L,,Mn, Nn}, satisfying Eqn 2.18. Presumably, any arbitrary solution of the diffrac- 
tion equation can be expressed as a linear combination of these vector functions. 
However, the existence of a generalized Fourier series expansion supposes that the set 
{L,,,M,,N,,} forms a complete set. 

Proving the completeness of the set {L,,M,,N,} is a two step process. First, 
one has to prove the existence of a complete set of functions for the diffraction equa- 
tion. Next, one has to show that the {L,,M,,N,} set can indeed span the solution 
space of the diffraction equation, or equivalently, the vector Helmholtz equations. If 
the labeling index n is not countable, then an arbitrary solution could be expressed 
as a generalized Fourier Integral of the basis functions with respect to the labeling 


parameter. 


2.4 Completeness of Vector Eigenfunctions 


Our goal is to be able to represent the solution of electromagnetic scattering in a 
series of vector eigenfunctions of the diffraction operator. The fundamental question 
we have to address is whether or not such a set is complete (in a strictly mathematical 
sense). In other words, whether the set of functions under consideration forms a basis. 
It can be shown that the solution space of the diffraction operator cannot be spanned 
by any finite set of functions. Assume momentarily that we have obtained only the 
L,, functions from the scalar solution, and that we have no knowledge about the 
existence of the M,, or N, functions. So we have a countably infinite set of functions 
that satisfy the diffraction operator. But it is easy to prove that such a set does 
not form a basis. In other words, there are elements in the solution space that are 
independent of the L, functions. For example, the M,, functions defined as L, x a are 
clearly orthogonal to the L, functions. Naturally, the same concerns are valid even 
with the knowledge of the larger set containing all the three types of functions. We 
do need to address the question of whether or not the set {L,,M,,N,} is a complete 
set. 

It is a well established fact that the set of all plane wave solutions form a com- 


plete set in the solution space of the vector Helmholtz operator. Physically this 
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implies that any scattered wave can be constructed by linear addition/superposition 
of plane waves. The diffraction operator is somewhat more general than the vector 
Helmholtz operator in the sense that it is satisfied by ‘generalized plane waves’ which 
have non-zero divergence. The vector Helmholtz equation, which follows directly 
from Maxwell’s Equations in a homogenous and isotropic medium is satisfied only by 
the zero-divergence solutions. Within a non-isotropic medium in which momentum 
transfer can occur, such as in a crystal, the E need not be perpendicular to k, and 
the zero-divergence solutions cannot describe such a wave. Although we are not di- 
rectly concerned with non-isotropic media, for the sake of arguments of mathematical 
completeness, we shall work with the solution space of the more general diffraction 
equation. 

Although the set of all generalized plane waves span the solution space, there is 
a fundamental difficulty with such a set. It arises from the fact that such a set is 
uncountable. In other words, the label(s) to identify the individual elements of the 
set are in this case, continuous variables (being the value of the propagation vector, 
and its direction cosines, and the direction cosines of the electric field). Such sets are 
not easily amenable for construction of general scattered wave solutions. It is here 
that a countable basis comes to our rescue. 

To begin, let us identify a few properties of the diffraction operator and its corre- 


sponding solution space. 


1. The space is linear. That is addition of two solutions and multiplication of a 
given solution by an arbitrary complex number is also a solution. This follows 
directly from the linearity of the diffraction operator. In particular, 0.f = 0 for 


all f belonging to this space. 


2. The space has an inner product and a metric. Mathematically this implies that 
with any pair of elements a and b, we can associate a complex number called 
the inner product (a|b) satisfying the following rules: 

e (calb) = c(a|b) where c is an arbitrary complex number. 
e 


d+ al|b) = (d|b) + (alb) for any d belonging to the set. 


e (alb) = (bla)* where * represents the complex conjugate. 


( 
( 
( 
( 


e (ala) > 0 for a 0 and (ala) = 0 only for a=0. 
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3. The space is complete since there exists the plane wave set which we know 
is a complete set. Mathematically it means that every Cauchy sequence is 


convergent. 


With these properties satisfied, the set of all solutions of the diffraction operator form 
a Hilbert Space. The inner product can be defined as will be seen when considering 
the solutions in a particular coordinate system, such as the spherical polar coordinate 
system (Chap 3). 

It can be shown that if {e;} is an orthonormal set in a Hilbert Space H, and if x is 
any vector in H, then the set S = {e; : (x|e;) # 0} is either empty or countable [37]. 
The importance of this theorem is that this guarantees countability of the basis if it 
exists. It can also be shown that every non-zero Hilbert space contains a complete 
orthonormal set. This follows from fundamental axioms of set theory embodied in 
Zorn’s Lemma [17]. The important fact to consider at this stage is that the two 
statements mentioned above concerning Hilbert spaces tell us immediately that there 


exists a countable basis. 


2.4.1 Completeness of the L, M, N set 


We have already seen that we can arrive at a set {Ln,M,, Nn}, of vector eigenfunc- 
tions that are mutually non-coplanar. It is clear that Vw, and Vw, are different 
functions in general. Similarly, it can be argued that the set of functions is such 
that at any given point in space, they are not all pointing in the same direction. In 
mathematical language, the set is linearly independent. The set also has a countable 
infinite number of elements in it. The scalar functions {v,,} are continuous with at 
least continuous partial derivatives up to second order. This implies that the derived 
vector functions are continous as well. 

Assume that we are given an arbitrary solution to the diffraction equation, x. By 
saying that we are given the solution, we mean that its value has been specified at a 
given set of points, S = {x1,X2,...}. For the moment let us assume that this set S$ 
is a countable set, as we have indicated by subscripting with natural numbers. For 
example, if the value has been specified at all coordinates (x,y,z) where x,y,z € 
{p/q: p,q € I}, then the set can be shown to be countable. Thus the specified set of 
points can be put in one-to-one correspondence with an index set, such as the natural 


numbers. 
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Now if we pick the first specified point where the function is defined, we obtain 
a vector of a certain magnitude which points along some specified direction. We can 
immediately pick the triplet of solutions {L1,Mj,N,}, and by virtue of their linear 
independence, we can choose suitable coefficients so that the sum a,L, + 6,M, + 
c.N, = X1, where the coefficients a1, 6,,c; are complex numbers. Since the left- 
hand-side is a linear combination of solutions to the diffraction equation, so we have 
a valid solution that is equal to the specified solution x at one point. In general 
the diffraction operator will propagate the linear combination a,L, + 6;M, + «Ni 
so that it will be different from the specified solution at other points (if it does 
not, then of course we have obtained the desired expansion, and we can stop the 
process here). So assume that the linear combination just obtained deviates from the 
solution at point 2, x2. Now we can pick our second triplet {L2, M2, No}, and obtain a 
‘correction’ to the original solution so that it matches at both the points. Essentially, 
we are solving a system of six equations in six unknowns to satisfy the match at 
the two specified points. Now it is easy to see that we can continue this process to 
‘match’ the specified function in an infinite series of the set of vector eigenfunctions 
to any arbitrary precision. By virtue of continuity of these functions, their linear 
combinations are also continuous for any finite number of terms. The difference 
between the arbitrary solution and the series just obtained can in principle be made 
as small as we wish. In other words the sequence of approximations will eventually 
converge to the desired solution (it forms a Cauchy sequence). As a corollary, we 
immediately note that the set of expansion coefficients so obtained need not be unique. 
It does depend upon the order in which we picked them to satisfy the conditions for 
a match. 

The validity of the statement of convergence of arbitrary Cauchy sequences follow 
from general considerations of a more restricted space of square-integrable functions, 
L?, with a semimetric. We know that the scattered solutions have to satisfy the 
radiation conditions. In other words, they have to vanish at infinity in a square- 
integrable sense. This follows from the finite energy content in any scattered wave 
from finite objects. The L, M and N functions satisfy such conditions. It can be 
shown (The Reisz-Fischer Theorem) that every Cauchy sequence in the semimetric 
space L* converges to a function in L?, or that it is complete [2]. 

To conclude this section, we note that the existence theorem concerned itself with 
an orthonormal set. The set of functions {L,,M,,N,} are not entirely orthogonal. 


However, they are a linearly independent set. This allows us to invoke the process 
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known as the Gram-Schmidt Orthonormalization procedure to obtain a complete or- 
thonormal set from the given set of linearly independent vectors [37]. In practice, as 
we shall see later, it is not always required to carry out the process of orthonormal- 
ization to begin with. 


2.4.2 Zero Divergence Solutions 


Consider a solution F whose divergence is zero. Let us find an expansion of F in 


terms of the basis {L,,M,, Nn}, so that 


F = SS {anM,, + bpNn + ¢nLn}. (2.26) 


Taking the divergence of both sides of the above equation, we find 


V-F So {anV > Mn + onV > Na tenV -La} 


\ {av La}: (2.27) 


nr 


or, 0 


Since this has to hold true at all points, we conclude that all the c, must be zero. In 
other words, a zero-divergence solution can be expressed only in terms of the M and 


N functions. 


2.5 Boundary Conditions for Time-Varying Fields 


Consider an interface between two media | and 2, with surface charge density o(r, t) 
and a current densities Ji(r,t) in medium 1 and J2(r,t) in medium 2. To avoid 
confusion with notation, conductance will be represented by g(r), since o(r) has been 
used to represent the surface charge density. Assume the volume charge densities 
to be pi(r,t) and p2(r,t) in the respective media separated by the surface. Let us 
consider a pillbox as indicated in the Figure 2.1. It has cross-sectional area AS and 
extends 6z, and 6z2 respectively into medium | and 2. The unit normal directed into 


medium 1 is denoted by n. From the continuity equation for charge: 


V -J(r,t) = a (2.28) 


Integrating both sides of Equation 2.28 over the volume of the pillbox, we obtain: 


[Vv -Se,t@r = ae | / p(r,t)d°r = eeu. (2.29) 
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Figure 2.1: Pillbox at an inteface between two media. 6z, and 6z2 are 
made to approach zero so that the only contribution to the suface integral 
occurs from the area AS only. 


where Ag(t) is the total charge within the elemental pillbox and is given by : 
Aq(t) = pi(r, t)AS é21 + po(r, (AS éz2 + o(r, t)AS. (2.30) 


If we let the pillbox thickness shrink, so that 6z1 and 6z2 go to zero in Equation 2.30, 


and also employ the divergence theorem, we obtain : 


fs-as=—5 [oe.nas (2.31) 
[ie I(r, 1)-w) ds = o(r,t)ds (2.32) 


which, in the limit of dS approaching zero, leads to 


do(r,t) 


(Ji(r,t) — Jo(r,t))-n = — i 


(2.33) 


By similar reasoning, since 


V+ Dr, t) = (r,t) (2.34) 


we can write the difference in the normal component of D at an interface 


(D,(r, 1) — Da(r, t)) - a = o(r, 2). (2.35) 
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For sinusoidal time dependence of the excitation field, the surface charge density 
o(r,t) = oo(r)e"’. Also, J(r,t) = g(r,w)E(r)e~’. So from Equation 2.33 we can 
write 

gi(r, w) Ey (r,t) - a — go(r,w)Eo(r, t) -@ = two(r,t). (2.36) 


Using Equation 2.35 this becomes 
gi(r, w)E, (r,t) -n — g2(r,w)E,(r,t)-n = iw(Dy(r,t)-n—Do(r,t)-n) (2.37) 
Using the constitutive relation D = cE we obtain 
(a+ “$1 B= (at “$2 BE, - A. (2.38) 
w w 


So, finally we obtain the time-varying boundary condition, the normal component of 


(complex) D is continuous 








Ey (r, t) n= €,E2(r,t) “ni. (2.39) 








Here € is the complex dielectric function of the corresponding medium. Thus we find 
that Equation 2.39 assumes the same form as the equivalent static case by using the 
complex dielectric functions. It is the complex dielectric function that incorpo- 
rates the losses and phase delays associated with time-varying fields in a medium. 
Consider a closed loop across an interface, as indicated in the Figure 2.2. By 


integrating both sides of the equation 








vine jee (2.40) 

ot 

over the surface enclosed by the loop, we obtain 
7 OB(r, t) 

wv x E(r,t)-dS = -f a es (2.41) 

0 
pE(r,t)-d = -— [Bas (2.42) 
¢ B(r,t) dl = 0 (2.43) 


since the surface integral will vanish in the limit of shrinking the height of the loop. 
From this it immediately follows that the tangential component of the Electric field 
(even in the time varying case) must be continuous across an interface between two 


media. Thus 








Ei(r,t)-¢ = Ea(r,t)- | (2.44) 
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Figure 2.2: When the rectangular loop in made to shrink in height, the 
contribution to the line-integral will be dominated by only the lengths 
parallel to the interface. 


If we exclude ferromagnetic materials, as well as materials with magnetic hystere- 
sis, then since we do not have magnetic monopoles or the equivalent of an electric 
current for magnetic fields, the boundary conditions for the magnetic fields will have 
the same form as the corresponding static boundary conditions. Specifically, the tan- 
gential component of H, the magnetic field, and the normal component of B, the 


magnetic induction must be continuous across a magnetic interface: 








H, (r,t) -t = Ho(r,t)- ¢ (2.45) 








and 





4H, (r, t) ‘n= fHo(r, t) n (2.46) 











where py and pg are the permeabilities of the respective media. 

Since there are two independent tangential directions and one normal direction 
for each point on a surface, these boundary conditions provide us with a total of six 
equations for each interface between two media. However, these equations may not 
be all independent of each other. It can be shown [22] that for a plane interface 
and an arbitrarily polarized plane wave, the equations on continuity of the tangential 
components of the E and H fields are the only independent equations. For non- 
planar surfaces, that are sufficiently smooth, the same arguments should hold locally 


at any point on the surface where the curvature can be considered negligible and the 


19 


arguments for a plane should hold good. Since this holds true for an arbitrary plane 
wave, it should hold true for any electromagnetic wave such as non-planar or even 
non-periodic waves. Non-smooth structures such as an edge can be be approximated 
by a surface with finite curvature. 

In general, obtaining the solution of vector fields imply that we obtain the solution 
of three mutually non-coplanar fields. On the boundary surfaces, we can construct 
three functions from the solution: Two of them representing the tangential compo- 
nents of the fields on the boundary and one representing the normal component on 
the boundary. The continuity of the three components across the boundary leads 
to unique determination of a vector field (with three components). For the electro- 
magnetic field, the E and H are coupled. So, although we can determine each field 
independently with three independent equations each on the boundary surfaces, the 
coupling of the fields tell us at once that not all of them are going to be indepen- 
dent. In summary, one can obtain four independent equations from the boundary 
conditions for each interface between two media for special geometries. For arbitrar- 
ily shaped boundary surfaces we need to consider all the six equations and solve an 
over-determined system of equations. 

Finally, one has to satisfy the radiation conditions which requires the scattered 
solution to decay inversely as the distance in the far field, and to assume the form of 


a spherical wave. 
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Chapter 3 


Plane Wave Expansion in Basis Functions 


3.1 Introduction 


The success of the method of (vector) basis function expansion relies on our abil- 
ity to express, conveniently, the incident radiation in terms of these basis functions. 
Although the completeness of this basis set guarantees the existence of unique coef- 
ficients for the expansion of an incident plane wave, in practice this involves some 
non-trivial algebra. It is our goal in this chapter to obtain explicit expansion coeffi- 
cients for a general plane wave in these basis eigenfunctions. Linearity of Maxwell’s 
equations imply that any excitation, non-planar or even non-periodic, could be solved 
in principle, from the known solution to plane wave excitations by use of Fourier space 
analysis. In principle, the set of basis functions obtained in any general orthogonal 
coordinate system could be considered a valid set of eigenfunctions in which the so- 
lution could be expressed. However, strictly algebriac and numerical considerations 
are in favor of the use of the functions pertaining to the spherical-polar coordinate 
system, for systems with azimuthal symmetry. The extensive analytical foundation 
readily available for the spherical Bessel/Hankel functions as well as that for the 
Associated Legendre functions allow us to obtain the expansion coefficients in closed 


form expressions. 


3.2. The L,M,N basis 


The solutions to the scalar Helmholtz equation in spherical polar coordinates are 


functions of the form (Appendix A): 
Pmnl?, 9,0) ~ 2n(kr)P™ (cos 0)? (3.1) 


where z,(kr) represents either the spherical Bessel functions j,(kr) or the spherical 
Hankel functions of the first kind h) (kr). The spherical Bessel functions are regular 


at the origin, whereas the spherical Hankel functions diverge at/near the origin. So 


Zi 


a region including the origin can only feature the spherical bessel functions in its 
expression for the field. A region not including the origin can have contribution 
from either of these functions. The labelling indices are n € {0,1,2,...} and m € 
{0,£1,+2,...,tn}. The P7™(cos@) are the Associated Legendre Functions. 

We have seen earlier that we can obtain a vector function L = Vz which is a 
solution to the diffraction equation if x is a solution to the scalar Helmholtz equation. 
Thus we can define a set of functions L,,,, in the spherical polar coordinate system : 


Linn = r : : 
“Op oe oe “oF sin O Og 








(3.2) 


Since the solution ~,,,, are available in separated form, the partial derivatives reduce 
to total derivatives since they act on functions of a single variable only. Thus we 


obtain the L,,,,, functions as 














dz,(kr) in(kr) dP™(cos@) , 
Lies = mn = k p™ 6 amd : n amd 
Vw d(kry ° (cos #)e'™*e, + fe 6 een 
_ 2n(kr) P™(cos 8) ae 
+em ia sind Og ys (3.3) 











We have already seen that the M,,,, functions can be obtained by defining M = 
V x aw, where a is an arbitrary constant vector. However, for the spherical polar 


coordinate system, we can obtain a set of M,,,, functions by defining 
M = V x (u(r)we,), (3.4) 
where u(r) is an unknown function that is to be determined. Thus we have 


e, reg rsindeg 





= a a a 
M r2sin@| or 28 a6 
uy 0 0 





_ 1 auy), 1 (Aur) 
—  rsind ad re = ( oe es Cy 


Now we need to impose the conditions that M has to satisfy in order to obtain the 


constraining relations for the function u(r). Assuming continuity of the functions and 
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their derivatives, we have 


: = ae Oia : oO (urb) 
MOND “Fe sat i d +5 O¢ 


+5 (22) 


1 fou) ue) _ 
r2sin@ | 000¢ 0¢00 f 





(3.6) 





Thus the divergence is zero for the M functions chosen according to Equation 3.4. 


We require the M functions to satisfy the diffraction equation 
V(V-M)-VxVxM+k’M =0. 


In order to obtain the necessary conditions to be satisfied, let us evaluate the V x V x M 


term. We have 








e, reg rsin deg 
VOM ee 2 
72 gin | oF ee eo 

1_ A(uy) + Balu) 

end Oe sin 0-5 


(3.7) 


_ 1 Oud) —._, 0*(urd) 1 0?(urp) 
sind | cont oe ane 06? sin@ Od? he, 


_ 0° (uv) (ur) 
+ {rind ar6 bert Drdd hey]. 

















Therefore, 
e, reg rsineg 
a a a 
; or 36 36 
VxVxM = and 
r? sin 
r2 a ( £08 gw) 
+ Ho (uy) a (up) —- ? (urh) 
sin O53 aro6 drab 








1_ 3 (up) ) 
sind d¢2 
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The e, component of V x V x M is 


1 [eed - Be 
r2sin@ |000r0¢ O¢0r00| 





(Vx VxM), = 


From Equation 3.5 we have the radial component M, equal to zero. Thus the radial 
component of the diffraction equation is satisfied by the current choice of M. 


The eg component of V x V x M is given by 
1 1 
(Vx VxM)s = g | 3 (sino?) 








rsinO?0¢| r?sin@ | 06 00 
1 Pup) (ue) 
r?sin?@ 0¢? Or? | 


So we can write 


(-V x VxM)-e9+k?M-e, = 











1 8 [A (ur) 1 @ (sino) 1 @ (ued) 


2 
rsinf0¢ | Or? r2 sin 6 06 Oo r?sin?@ 0¢? Fe ww} 


In order for the diffraction equation to be satisfied, the quantity within the square 


brackets must be independent of ¢. By following the same procedure, we obtain 


(Vx Vx M)-e,—k?M- eg = 








10 |? (udb) 1 @ f/f .-,d(uy) 1 0? (uw) 

6 i : 

r00| Or? rsind00 \ 00 r2sin?@ 0¢? a OY 

As we can see, the expressions within the square brackets for both the @ and ¢ 
equations above are one and the same. Thus the diffraction equation will be satisfied 
if the quantity within the square brackets above is some function of r alone, say f(r). 


In that case the right-hand sides will vanish. In particular, we can choose f(r) = 0, 


so that the condition on u(r) is given by 


O*(urp) 1 Of. ,A(up) 1 Pub) 1. 
OF r? sin 8 00 pine 06 + r2sin26 Od? + k*uy = 0. 











In addition, if we choose u(r) = r, we obtain 


(rv) 1 a (sino) 1 aa 








kr = 0. 
Or? rsin@ 00 Oo r sin? 6 0¢? ae 
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The left-hand side of the above equation is equal to r- (V7 + k?w). So for the above 
choice for f(r) and u(r), we have the diffraction equation satisfied by all the three 
components of M, when 2) satisfies the scalar Helmholtz equation. Thus we have : 
Mina: = Vo (rein) 
= Wmn(V x re,) + Vtban X Ter 
So Waa eT ( Since V x r = 0) 
= Lay xr 


Therefore, we obtain the desired M.,,,,, functions that have zero divergence : 





P™ (cos 8) 


sin 8 


Minn = tmz, (kr) ete, — Zn(kr) 


dP™(cos@) ,., 
—9¢ ey. (3.8) 











Having obtained M,,,,, we can easily obtain the other set of zero divergence func- 


tions, the N,n», since 


1 
Nk = i x Minn 


k 
e, reg rsin Jeg 
= a nee a a 
kr? sin @ | 9 oF oP 


0 +r (imzn( kr) FECA) cine) r sin @ (—zn( br) PEG) cis) 





1 _ 0? P™(cos 6) OP™ (cos 6) » P™ (cos 8) 
= Gagan r=attr) sin 8 762 cos 8 70 +m ae a e, 


is (rt) AP (cos 4) Sra(tr)) Pr"(cos Me, See 
r 





nO” 
a Or sin 8 


06 


Since P(cos @) is the solution to the 6-part of the separated scalar wave equation, 


eg + am ( 


therefore 
i eke cos 6 d.P™ 1 m? p™ =0 
sae ae ae 








Equivalently, 
0? P™ (cos 9) OP™(cos 8) » P™ (cos 8) 
06? ey oe ree sin 8 





sin 6 = n(n + 1)P"(cos 8) sin 4. 
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This simplifies the e, term and we obtain, remembering to change the partial deriva- 


tives to their corresponding total derivatives : 














me ; lvd(re, (enya Oy. sa. 
ie. ne 4 1/2 (EY) pr(eos Oe’ %e, a = Wrenk r) (cos ) ei te, 
_ 1 d(rz,(kr) P™ (cos 4) ae 
+im 7s a sind ee ey: (3.9) 











It is of some interest to note that since the expressions for Linn, Mmn and Nin 


have the z and m occuring as a product, we can easily derive the following useful 











relations: 
(nm)! 
Minn = (-1) (mn tm) vie (3.11) 
m(r—myl, 


3.3. Orthogonal Properties of L, M and N Functions 


Having obtained explicit expressions for the vector functions L, M and N, it is desir- 
able to examine their orthogonal properties. Denoting the complex conjugate of L,,, 


as L* 


x ny we have 


dz,(kr) dz* 
d(kr)  d( 





mint 


(kr) 
Al. 
(3.13) 
dP™(cos@) dP™ (cos8) — mm’! ; 2n( kr )z*,(kr) 
7 Weak 6) P” g) | — +“ |. 
+( do a ag ee (kr)? 


| Bere Ba = fieimee-imis [Petco 0) P” (cos 8) 





On carrying out the ¢-integration, within the limits of ¢ € [0,27], since 
20 : a, 
| dei e-i™'6 = nF amt (3.14) 
0 


We obtain 
dz,(kr) dz*,(kr) 
d(kr)  d(kr) 





27 
i db Limn- Ltn = 20k? [r(cos 0) F(oos) (3.15) 
0 


26 





dP™ (cos 8) dP? (cos 8) m? in(kr)z*, (kr) 
= n pm 6) P” g) | —~_~4 + |. 
( dO age age ge Se a (08) es 


Using the orthogonality relations for the Associated Legendre Functions, 


1 2 (n+m)!)! 
dxP™(x)P™(x) = Rie i 
I aE eee) 2n+1(n—m)! 1) 





and the following relation that we derive in Appendix C, 


To apd dP" (cos 8) dP™ (cos 6) Ware 28s a 
‘ dO sin 6 ( 6 6 + a2 gin (cos 8) P'")(cos 8) 
2n(n+1)(n+m)!)! 
— Oat: wl 
2n+1 (n—m)! (317) 








we obtain, on carrying out the 6-integration, 


27 wT 
| dé if d6sin@ Lian Legs (3.18) 
0 0 


2n(n +1)(n+m)!)! [eta n(n +1) 
2n+1 (n—m)! d(kr) (kr)? 


where the square denotes the square of the absolute value. Using the recurrence 





= 2k’ iain | Online 


relations on the spherical Bessel/Hankel functions, with p = kr : 


dzn(p) 1 
dp ~ On+1 





{n2na(p) — (n+ 1)2n4(0)}, (3.19) 


and 





= = = [ Ven-1(p) + 2nt1(0)} (3.20) 


we obtain, keeping in mind that the squares of the Bessel/Hankel functions denote 


the squares of their absolute values in these expressions: 





20 wT 
j do | Ge ati 
0 0 


= re Ee ctalke) + (mt Del] fob 82D 














Similarly, the orthogonality relation between the M functions can be obtained: 


Zt 





apn 4 1 ! 
/ dg | M6 Gee Aa NA gl a ged, 163) 
0 0 2n+1 (n-—m!)! 














Using the following recurrence relation on the spherical Bessel/ Hankel functions: 


1. 
Er dp tent r)) 





= aa {(n + Lanai (kr) = Nin4i(kr)} ) (3.23) 


we obtain the orthogonal relations for the N functions : 





20 Tv 
7 dd | dO sin @ Naan « N* mgs 
0 0 


= A ng 1)h b+ mcte lb] Soba (828 














Now we need to examine the cross terms; orthogonality relations between the L, 
M, N functions. The obvious ¢-integrations tell us immediately that functions with 
different m’s will be orthogonal to each other. So we need to examine the behavior 
of these cross terms for differing n’s only. We can write 
gn(kr)z*,(kr) 1 d 

‘f sin @ dé 
For m = 0, the expression on the right-hand side vanishes. For m # 0, the 6- 





em (P (cos 6) P77" (cos 6)3.25) 


20 
i dé Lng Mio = —Daimk 
0 


integration vanishes since P(x) = (1—2?)d™P,,/dx™ goes to zero at the limits when 


x? = cos*§ = 1. Thus we have the orthogonality relation (for all n,m): 





20 wT 
fi dé if d0sin@ Lmn-M%y =0  (Vn,m,n',m! eT) (3.26) 
0 0 











By reasoning in exactly the same way we can show that the M- N* term vanishes 


for all values of n,m : 





20 ™ 
i do | désin8@ Minn: N* 1. = 0 (Vn,m,n',m' € I) (3.27) 
0 0 
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Now we examine the cross term L- N: 


20 T 
y dd i; dO sin @ Ln Nein 
0 0 
Boil RE) den (kT: 
kr d(kr) 





= 2k fe +1) i dé sin @ P’" (cos 6) P-"" (cos 8) 
0 


rh a 2 
er) tf (re (hr)) 


OF aca dP™ (cos @) dP"}(cos 6) Hs wea. os 
i d@ sin 0 ( 10 10 + ant gin (cos 8) P*(cos 8) | > drm 








By virtue of Equation 3.16 and Equation 3.17 we can see that the 6-integrals will 
vanish for n 4 n’. The only isolated case when the integral may not vanish is when 
n = n', when we can use the recurrence relations on the spherical Bessel/Hankel 


functions to obtain: 





20 wT 
i dd | d6sin 8 Linn « N*y.) 
0 0 


= a ms . * = len_a(br) — zha1(br)] Sant Sram (3.28) 














Since the right-hand side of the above equation is a real number, if we obtain 
the complex conjugate of either side of the equation, we immediately come to the 


conclusion that 





20 wT 20 wT 
f dé | dO sin 8 EaciNige | dé | d0sin0 Ninn Lt, (3.29) 
0 0 0 i) 











Thus, the L,M,N basis functions are not completely orthogonal to each other. 
The overlap between the L and N functions for the same n requires us to consider 
both the matrix elements even for cases when we are trying to expand a plane wave 
that has zero divergence. However, we can easily observe that the {Linn, Minn, Nn} 
functions are mutually non-coplanar. Since they also form a complete set of functions 
in the Hilbert space of all solutions to the diffraction equation, in principle one can use 
the Gram-Schmidt Orthogonalization procedure to derive an orthonormal set from 


the given set of functions. However, as we will see shortly, such a theoretical process 
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is not required for most applications. We can use the L,M,N basis functions directly 


for many types of computations. 


3.4 Plane Wave Expansion in L, M, N Basis 


The completeness of the set of functions L, M and N assures us of the existence of a 
valid expansion series for an arbitrary plane electromagnetic wave which satisfies the 
diffraction equation, in this basis set. Once a valid set of coefficients are obtained, 
we can offer an operational proof of the completeness theorem. If any plane wave can 
be represented in a convergent series of these functions, then it immediately implies 
that the set must be complete. This follows from the fact that any scattered solution 
to the vector Helmholtz equation can be constructed from a plane wave basis. 


To begin our discussion, let us consider an arbitrary plane wave: 
F = Ee‘**. (3.30) 


Here E = E,e, + E,e, + Eze, is the (oscillating) electric field of the plane wave, 
k = (k,a,) is the propagation vector and r = (r,6,¢) is the position coordinate. 
To keep our derivation sufficiently general, we shall allow the orientation of E and 
k to be along arbitrary directions. In reality, for a plane electromagnetic wave in 
a homogenous and isotropic medium, the E and k must be orthogonal, and so the 
divergence of F must vanish. However at this stage we do not need to impose such 
constraints. We shall also allow the E,, E,,£, and k to be complex numbers. This 
represents any general elliptically polarized wave propagating in a medium that could 
be attenuating /absorbing or even amplifying. 


The completeness of the set of functions validate the following representation: 


ee = la? Mie POL Ned +c be} (3.31) 
eye = y) fat Mie be Nea + ba} (3.32) 
ee SS fae Mees be Nag Pc bas} (3.33) 


Therefore, we can write our general plane wave expansion as: 


Ee'** = an {(Ee@een =F Ey Gan a Eas) Minn 
+ (6%, + £0, + £81) Nn 


yrmn 


+ (yet, + EyGn + ExGin) Linn (3.34) 
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If E is represented by (F, 6;, ¢;), in the spherical polar coordinate system, then 


Ei, = Esin @; cos ¢; E, = Esin 6;sin ¢; E, = Ecos9§;. (3.35) 
Therefore, we need to obtain the nine coefficients {a7,,,, a.) ---; an} Which are func- 


tions of {n,m,a, 3}. 


3.4.1 Deriving a7,,, 6%, and c,,, 


mn? “mn 


It can be shown [43] that the exponential part of the expression for a plane wave can 


be expanded in the following series: 





eikr 


Me 


= (28-4 Libr) {- Bo 


0 


s 


Silos a)e"® P!(cos jeri | (3.36) 


Let us denote e.e*™ = E., where c denotes x, y or z. Let us introduce the notation: 


0° Tw 20 
| ar [ do sind f dé E,-X*,, = [Be Xin = (E.|Xinn) (3.37) 
0 0 0 


where X,,, represents Linn, Minn ot Nin. We can view (E.|Min,) as the definition 
of the inner product in the Hilbert space of all solutions to the diffraction equation. 


Eiko Ce.) ising. the 


From Equation 3.33, we can calculate the coefficients {a2,,,, 2.15 Gan 


orthogonality relations already obtained. Obtaining the inner product of both sides of 


the equation with respect to Linn, Mmn and Nmn respectively, we obtain the following 


relations: 
(Ez|Minn) = Gnn{Mmn|Mmn), (3.38) 
(Bz[Linn) = Bigg NonnlEmn) + ey (Lonn Line) (3.39) 


which can be readily solved to yield: 





z (Ez|Minn) 


mn = (Monn| Minn) (3.41) 
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2 (Ez |Lrnn) Lmn | Nin) — (Ez |Nimn) (Linn [Linn) 
Ban = TS A CT ee eT ee PT (3.42) 











(Ez|Nmn)(Nmn|Lmn) — (Ez|Lmn) (Nmn|Nmn) 
PET CIN ca Mien Vien Newnes ens) (Lane Lien) (3.43) 














By substituting the label z with x or y in Eqn 3.41, Eqn 3.42 and Eqn 3.43, we 
immediately obtain the corresponding formula for the remaining coefficients in terms 
of their inner products with the basis set functions. 

The integrals represented by the various inner products in the expressions for the 
coefficients require a good deal of messy algebra and careful analysis for its evaluation. 
The integrands involve products of spherical Bessel functions, Associated Legendre 
functions as well as their derivatives. The non-triviality of these integrations is am- 
plified by the fact that one has to evaluate a double infinite summation (over s and 
1) to arrive at a closed form expression for these integrals. The inner products be- 
tween the basis functions are obtained from the corresponding orthogonal relations 
Equations (3.22, 3.24, 3.21, 3.29). By carrying out the r-integration on both sides of 
the equation we can show (worked out in the Appendix D): 

2n°*k = (n+m)! (4n?+4n—-1) 


Linn Dran) = (2n + 1)? (n — m)! (2n — 1)(2n + 3) nn 





Maan) = Eee os 





Qn? (rn +m)! n(n + 1)(4n? + 4n + 3) 
k(2n + 1)? (n —m!)! (2n — 1)(2n + 3) 





(Ninn [Nmn) (3.46) 


8x7 (n +m)! n(n + 1) 
(2n + 1)? (n — m)! (2n — 1)(2n 4+ 3) 





(Linn |Ninn) (3.47) 


Since 


e, = cos fe, — sin eg, (3.48) 
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the integrals involving e, will contain only terms of the form P”, where n’ can be 


shown to be equal to +1 only. The integrals (E,|Minn), (Ez|Nmn) and (Ez|Lmn) 


therefore have the follwing form for m > 0: 





Qmr2artt 
E Minn = —__P™ imp A 
(E:IMo) = ZEE Pm(cos a (3.49) 
a9 2,n—-1 -impB 
(E,|Nmn) = 2 - (3.50) 


k(2n + 1) (Qn —1)(2n +3) * 
[n(n —m+1)(2n — 1) P™,(cos a) —(n+1)(n+m)(2n + 3)P™ (cos a)| 


n-1 


oa aaa ewe 
(Qn +1) (Qn — 1)(Qn +3) * 
[(n + m)(2n + 3)P™ (cosa) + (n —m+1)(2n — 1) P77,4(cos a)| 





(Ez|Lmn) (3.51) 


n-1 


On substituting the expressions for these integrals back into Equations 3.41, 3.42 
and 3.43, we finally obtain the functional form for the coefficients (for m > 0): 


























ann, m,Q, B) = gett m(2n at 1) (n = ME Spite aje i? (3.52) 
(m>0) n(n+1) (n+m)!)! 

b(n, m, a, 3) tl (n—m)! iis n(n —m + 1)P™ j(cos a) (3.53) 

= e€ . 
(m > 0) n(n+1)(n+m)! —(n+1)(n + m)P™ ,(cos a) 

(ny mya; B) 2-1 (n—m)! = 

a = <2 + 1)e7 0 ee 54 

fe) 7 aan n+l)e cos a P’ (cos a) (3.54) 











We observe that when a = 7/2, so that E and k are perpendicular to each 


z 


z , coefficients vanish, since the 


other as in a plane electromagnetic wave, all the c 
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divergence of such a field is zero. When a = 0 or g all the az,,, coefficients vanish 
since P™(+1) = 0 for m > 0, and a7,,, vanish for m = 0 because of the factor m in its 
expression. Similarly, the 67,,, vanish for a = 0 or z since P™(+1) = 0 for m > 0, and 
for m = 0 P,(+1)(—1)” for m = 0 so that the expression within the square brackets 
become for m = 0, n(n + 1).1 — (n+ 1).n = 0. So when a = 0 or x only the c?,,, 
coefficients can be non-zero. 

When m < 0, we can arrive at the corresponding coefficients by examining the 
changes that occur in the expressions for the individual inner product terms. As 
shown in the appendix, from considerations of parity of the Associated Legendre 


functions, we can arrive at the corresponding coefficients when m < 0: 

















a®_,, (n,m, Q, wind Sage (eI. 
7 . 0) 5) = (—1) cs e€ Ze = il Mm, a, B) (3.55) 
62 (m,m, a, af Sigg (aan s 
Ot ino Or savaene ene) 
= nee 0) 3) = (—1) € oe nl mM, a, B) (3.57) 











3.4.2 Deriving a%,,, 6%, and c,,, 


mn? “mn 


The process of finding the x-coefficients are similar to that of the z-coefficients that 
we have just obtained. However, the algebra becomes even more messy, because now 
we do not have the advantage of having the symmetry about the z-axis. Again, the 
detailed analysis is shown in Appendix D and we highlight only the key results from 


such calculations. We have 
e,, = sin 9 cos ge, + cos # cos Peg — sin Peg (3.58) 


So now, the integrals of the type (E,|L,,,), (Ec|Mimn) and (Ez|Nmn) will have terms 
with m+ 1 because of the extra sin ¢ and cos ¢ terms. The parity of all the parts of 
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the integrands considered together(the r,@ and ¢ parts) allows only the /! = m+1 


terms to survive. So form > 0: 





























(E,|Mmn) = geet (n+m)(n—m+ IP -*(cos ae en? (3.59) 

k(2n + 1) +P™t!(cos ae (m+1)p 

qeget | 
E.|Nmn) = : 

(Bo|Nmn) k (Qn+1)Qn—1)(Qn+3) ~ oo) 

(2n + 3)(n + 1)(n+ m)(n+m— 1)P™7"(cos a)eW"l 

—(2n + 3)(n + 1)P™+"(cos aj“ 
+(2n —1)n(n —m+2)(n-—m+t Ly Pra (cos a)eWl = 
—(2n — 1)n PA! (cos a)e~lmtt 
qeget 

cole aaa CTS CTEMEY ARTE) ie en 

(2n + 3)(n+ m)(n +m — 1)P™7"(cos a)je a ca 

—(2n + 3)P™* (cos a)e~lm41 
(2n —1)(n —m+ 2)(n-—m+t 1)P™5'(cos a)e7 -i(m— 
+(2n —1) PM (cos a)e~lmtt 
Using these expressions we obtain: 
x = ! 
az ,(m,m, a, B) _ nea (2n +1) (n-—m)! (3.62) 
(m > 0) 2n(n + 1)(n+m!)! 
(n+ m)(n—m+1)P™-!(cos a)e—tr—1)8 
+P™+1 (cos a)e~Hr41)8 
br (n,m, a, B) ah 1 (n —m)! : (3.63) 
(m > 0) 2n(n + 1)(n+m!)! 


(n+ 1)(n+m)(n+m— PP eer 
—(n + 1)PRY 
+n(rn—m+2)(n—m+t iD reare 
Sar 


COs & 


COs & 


COs & 


a i a a 


COs & 





SS 
! 
3 
| 
e 
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cz (nm, m, a, B) 5 wt (n —m)! . (3.64) 
(m > 0) 2k (n+m!)! 
(n+m)(n +m —1)P 7'(cos aje re 
—P™*1(cos ae +1)8 
—(n—m+2)(n—m+ 1) Pi, (eos ae Hr—1)8 
+ Pmt} (cos a)e~Hrt4)8 











For m < 0, the easiest way to obtain the coefficients (explained in Appendix D) 


are the following recipes: 














Cini, See : as = ke with ie cots (3.65) 
changed to e(™+)8, 

(3.66) 

pe 41)? (n +m)! is with eWtm™t1)6 ene (3.67) 
ye (n —m)! [changed to e(™#08, 

(3.68) 

oo = (-1)" (rn +m)! i with e707 =Ue nes (3.69) 
cae (n —m)! |changed to ei(™£08, 











by, and cY, 


mn “mn 


3.4.3 Deriving a’ 


The y-coefficients are obtained almost identically as compared with the x-coefficients. 
We realize that the only difference between the two sets arise primarily because of 
the interchanging of the cos ¢ and the sin ¢ terms due to e, being replaced by e,. We 


have 


e, = sin @sin ge, + cos Osin ¢eg + cos Peg (3.70) 


As explained in detail in the appendix(D), the ¢-integration will cause a change in 
the sign of the m+ 1 terms while the m — 1 terms will continue to have the same 
sign. Also, an extra factor of —2 is generated because sin d = (1/27)(e’® — e~"®). Thus 


following this general recipe we obtain the y-coeffcients: 
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Brn(%, M08) _— oy (2R+1) (n= m)! (3.71) 


(m > 0) 2n(n + 1)(n+m!)! 


(n+ m)(n—m+1)P™-!(cos een ee 
—P™+1 (cos a)e~™+1)8 











by (n,m, a, 3) i} (n —m)! 














ae 372 

(m>0) " On(n +1) (n +m)!” 2) 
(n+ 1)(n+m)(n+m— 1)P™>"(cos ae ir—1)8 
+(n + 1) P41 (cos ae t+H)6 
t+n(n —m + 2)(rn —m + 1)P™5"(cos a)eWtm—1)6 
+n P™*1(cos ae tr +18 

Cmn(@™,a,8) _ a (n—m)! » (3.73) 

(m > 0) 7 2k(n+m)! ; 

(n +m)(n + m—1)P™71(cos ae“ 4-8 
+P™* (cos a)e~Hr41)8 
—(n —m+ 2)(n —m+ eae AEA ‘(cos aje—tm—-)8 
— Pmt} (cos ae t1)8 











Transformations similar to what was used for the x-coefficients for m < 0 are going 
to be valid for the y-coefficients as well. Thus the y-coefficients for negative m can be 


obtained as follows( note the extra (—1) factor as compared to the x-transformations): 
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a’ = (-1)™ (n +m)! bs with e—tmt1)6 oe oT 
— (n —m)! [changed to e(™#08, 
3.75 
Winn = (-1y"* (n +m)! i with e~#¢™+1)6 oan) 


(n—m)! changed to e(™+08, 





y ( pet (n+m)! ie with e—(m+tl)6 mao 
Cimn = 


(n—m)! changed to e(™+08, 











3.4.4 Numerical Convergence of Basis Function Expansion 


It may be pointed out that these coefficients for the expansion of an arbitrary plane 
wave in terms of the L, M, N basis are not necessarily unique. Alternate sets of 
expansion coefficients can be obtained by, say rearranging of the basis. This follows 
directly from the general theory of orthonormal basis in Hilbert spaces [37]. Since 
one can expand any plane wave in this basis, as discussed earlier, one can obtain an 
operational proof of the completeness of this basis. It is also important to realize 
that simply proving completeness of a set is not the last word in computations. 
One has to address the question of convergence as well. So numerical verifications 
are absolutely necessary to validate a certain basis set. To summarize the results, 
convergence of better than 1% is normally achieved by taking terms up to index n, 
where n ~ 1.4|kr|. This can provide us with the guideline on how many terms to 
include for a problem in which the geometry can be measured in units of A, the 
wavelength of the electromagnetic field in question. 

In Figures 3.1 and 3.2 we show the expansion of arbitrary electromagnetic waves 
in the L, M, N basis. We have specified an arbitrarily directed wave vector k 
with direction cosines a and /, with the electric field oriented along some arbitrary 
direction specified by the angles 6, and ¢,. The values of the exact expression for 
Ee‘*™ and the series expansion in the L, M, N basis are compared along some 
arbitrarily specified line directed along (@,¢). We observe that the convergence is 
very good for |kr| < 0.75n and this is fairly independent on the choice of the specified 
directions for E and k as well as the path along which the comparison calculations 


are done. In Figure 3.1 the wave vector is real where as in Figure 3.2 the wave 
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Figure 3.1: Demonstrating completeness of basis function expansion. The 
expansion of an arbitrarily chosen wave in the L, M and N functions. 
The solid lines show exact values. The broken lines are computed from a 
truncated series in the basis functions. The above calculations are done 
with 15 terms (n = 15) for kr € [0,12], 6 = 37 degrees, ¢ = 59 degrees, 
a = 123 degrees, 6 = 83 degrees, 9. = 49 degrees and ¢. = 21 degrees. 
Convergence is very good for |kr| < 0.75n. The horizontal axes are in units 


of |kr|. 
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vector is complex. Complex wave vectors correspond to an attenuating or amplifying 
medium. We therefore demonstrate that these functions can be used to expand an 
electromagnetic field in a general medium. Since the specified E and k are not 
orthogonal in general (as in the above two cases), it is therefore possible to expand 
waves that are more general than plain electromagnetic waves in a non-attenuating 
medium. Later we will also examine the possibility to obtain the solution of Maxwell 


equations in the presence of sources, using the L, M, N basis. 


3.5 Numerical Evaluation of Lyn, Minn and Ninn 


Numerical evaluation of Linn, Mimn and Nmn require us to obtain the values of func- 
tions of 8 such as 4 P™(cos ), alee) and functions of r such as £4 (rz,(kr)) and 
20) It is clear that one has to exercise caution in evaluating these functions close 
to @=0or x and r = 0. In this section we indicate numerically stable computational 
equivalents to these functions. 

From Equation C.11, by adding the third and the fourth equations we eliminate 


the «P(x) term, and so we obtain: 








*pr(x) =5 a (n= m+ 1(n+m)P™ (ax) — PMH (x)], (3.79) 
and since 
# pm (cos 6) = _/T—at P™(z) (3.80) 
do” ap te 
Therefore: 
do fi cits re 
apn (cos 8) = 5 [prtt(a2) —(n =m +1)(n+m)Pr"(a)] . (3.81) 











From the definition of the Associated Legendre Functions (m > 0): 








m d™ 
Pea) = (yn 2)F Pye) 
m am 
= (—1)2 sin” 9—P,(z), 3.82 
(—1)8 sin” 0-2 PCa) (3.82) 
where P,,() are the Legendre Polynomials and x = cos 6, we immediately obtain: 
P™ (cos 8) - cee ee: 
0 = (—1)? sin 0 mb n(cos 9). (3.83) 
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Figure 3.2: Demonstrating completeness of basis function expansion when 
the specified wave vector is complex. The solid lines show exact val- 
ues. The broken lines are computed from a truncated series in the basis 
functions. The above calculations are done with 15 terms (n = 15) for 
krmaz = 10+ 8, 6 = 35 degrees, ¢ = 42 degrees, a = 97 degrees, 3 = 82 
degrees, 6. = 21 degrees and ¢, = 79 degrees. Convergence is very good for 
|r| < 0.75n. The horizontal axes are in units of |kr|. 


4] 


When m # 1, clearly the right hand side is zero when 6 = 0 or 7. So, when m = 1, 
P}(cos 9) d 








Mcoe = (a), (3.84) 
Fron Leeeudie & cunterential equations 
ia) Pie Pape (3.85) 
ig de 
socio 
P,(2) a — 5 P,(o) e mn FD) (a) (3.86) 
= 5-P22)+ mnt (a) (Ushi Sos thodeie, oF Pla): 


When « = +1, then P?(x) — 0. Thus: 





Le 6 
TT cos ) 0 whenm#Zl (3.87) 
sin 
1 
= — when ¢ =1,m=1 (3.88) 
1 
= (ap Met) when «= —1,m=1 (3.89) 
1 
=. - where = lone = 1 (3.90) 
2 
1 
= (U5 when « = —1,m=-1 (3.91) 











The functions involving the spherical Bessel and spherical Hankel functions of the 


first kind for kr — 0, are best evaluated by using the following recursion relations: 
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ae = ot 1) [en-1(kir) + 2n4i(kr)| (3.92) 
Tenlkr) = k 1 k 3.93 
d(kr) Bal ir) = (2n +1) [nzn—1( r) = (n + zn41( r)] ( ; ) 
ld i 
aa (rz,(kr)) = naa) [(n + 1) zn-a(kr) — nzngi(kr)]. (3.94) 











Here, the z,(kr) represents either the spherical Bessel or the spherical Hankel 
functions of the first kind where n > 1. The argument kr could be complex in 
general. We do not require to evaluate these terms for n = 0 since the L, M and 
N functions are zero for n = 0. The absense of 1/kr factors on the right hand side 
allows straightforward evaluation of these functions as kr — 0. It is of interest to 
note that the field on the z-axis will be given only by the m = 0,+1 terms, due to 


the vanishing of the other 4 P™ (cos 6) and Fi (cos @) functions on the z-axis. 
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Chapter 4 


Shell Problem 


4.1 Introduction 


As an application of the method of solution by expansion in vector eigenfunctions, we 
consider a problem of current interest. We shall solve the electromagnetic problem 
of scattering from concentric spherical shells. There is a lot of interest to understand 
the electromagnetic response of nanoparticles that have a shell or coating of some 
other material. We shall solve the problem for a single shell. However, it would be 
fairly obvious that the method can be generalized to the problem of multiple shells. 
An underlying assumption of this method is that the macroscopic dielectric function 
of a particular material continues to be applicable in the nanometer regime as well. 
For clusters larger than a few hundred atoms and for shell thickness greater than a 


few monolayers, this is usually a valid approximation. 


4.2 The Model 


A sphere is centered on the origin, of radius R, of a material with dielectric function 
€4(w) ( Figure 4.1). Surrounding the sphere is a concentric shell of radius Ry of a 
material with dielectric function €2(w). The dielectric functions are complex functions 
of frequency of the electromagetic radiation in question. The sphere and shell are 
embedded in a medium, the dielectric function of which is specified as €3(w). Without 
loss of generality, we shall assume that we are dealing with non-magnetic materials, 
so that the relative permeabilities of the three region are approximately unity. 

A plane wave of unit intensity is incident along the z-axis, with its polarization 
along the x-direction. By symmetry, it doesn’t make any physical difference in the 
choice of the incident direction. However, the expansion coefficients for the plane 
wave in the {M,N} functions assume simple expressions when the wave is incident 
along the z-axis. In this coordinate sytem k is specified as (k,a = 0,3 = 7/2). Thus 


Gmn = ay, and by, = bf. Since the divergence of a plane wave is zero, we can 
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Figure 4.1: Shell-geometry with single shell. There are three spherically 
concentric regions indicated. The core has dielectric constant €,, the shell 
has dielectric constant €,; and the exterior has dielectric constant ¢3. In 
general there can be more number of shells. 


expand the incident field in the M and N functions alone. The scattered waves are 
the ‘forced’ solutions, and they will assume the same form as the incident waves at a 


given frequency, and so they too can be expanded in the M and N functions only: 
Now for a plane electromagnetic wave H = ay x E. Therefore: 


H 


1 
—S0 {aimnV X Minn + bmaV X Ninn} 
1 


Komn 


Using a = 0 and 8 = x/2, only the n > 1 and m = +1 coefficients turn out to be 


non-zero. We obtain for the expansion coefficients: 

2n+1 
2n(n + 1) 
n(n + 1)ain (4.4) 


x 
Gin 


a nel Oe) 2) 4" 


x 
Gin 


2n+1 
2n(n +1) 
By, = —n(nt+1)bip (4.6) 


i = Boe Laat 


Let us represent the core as region 1, the shell as region 2 and the exterior as 
region 3. Let us assume the following form for the solutions of the scattered E field 


in the three regions: 
B= 5 {oll MH + of NGL + oh Mah + Ge NC 
n=1 


The +1 indicates the two choices for the value of m. The superscripts identify the 
region and the type of function being used. For example, the A superscript in mM!" 
represents the M functions that are derived from the scalar solution that uses spherical 
Hankel functions of the first kind for their radial part. The 7 superscript represents 
functions that have spherical Bessel functions for their radial part. The numeric 
label identifies the region and so the value of & at the frequency of the incident 


electromagnetic wave. Specifically, we have (with m = +1): 


MMs = irnja( kyr) A) eintey — jn(kyr) AED) pind, (4.10) 

Ms = imma har) 2S) cimtey — j, (bar) AID) ind, (4.11) 

MOQ! = iml!(hyr) ALD cimtey — HY (byr) AICS cine, (4.12) 

M&)s = imnjal kar) 2) imtgy juke) ELE) cine, (4.13) 
P™(cos6) dP (cos6) 


M@&)h = imh (kgr) ee, — hY (kar) ee (4.14) 


sind " dé 6 


In( Kir) 
kyr 





Ns = n(n + 1) P™(cos6)e'”*e, 
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1 d(rjn(kir) [dP™(cosé) , F . P@(cos6) , é 
nr am nr am 4] 
5 ae ee me Se) oe 


Ns = n(n + 1) 


—In( k . 
JolBar) 27) Bee os0)e"*e, 





1 d(rjn(kar) { dP7" (cos) : P(cos0) , é 
el em 4] 
se dr dé cane sind ne (E18) 


AY (kar 
N@)A — n(n + 1) 2") pm (cost moe 
ar 














- arti (ha) [Pelee be apg cme) ren 

N®F = n(n +1) nt) pe (co s0)e™e, 
+ ~ Ae dolisr) (= oost) ci™bey + ima eos’) icos6) cntes| (4.18) 

NG = nin + 1) C9") po eosd inde, 
ee drift) («re Cor ee eel) cnt) bay) 


Since region | includes the origin, only the spherical Bessel function solutions are 
allowed. The spherical Hankel function solutions are singular at the origin. In region 
2, which does not include the origin, we can therefore have both the spherical Bessel 
and spherical Hankel function solutions. Region 3 does not include the origin. So we 
assume an expansion in the spherical Hankel function solutions only. The spherical 
Bessel function solutions are included in the expansion of the incident field in region 
3. Far away from the origin, where the scattered solution must become negligible, we 
must recover the plane wave field of the incident wave only. Thus this regularity at 
infinity restricts the scattered solution in region 3 to have only the spherical Hankel 
function solutions. The incident field itself has to be regular at the origin (since 
the choice for locating the scattering object at the origin was purely a matter of 
convenience) and so it is expressible only in terms of the spherical Bessel function 


solutions: 


=> {ai1, M Se ra Vein No} (4.20) 


AT 


The expression for the magnetic field intensity H in the various regions are: 


H, = fey fais, NEY + by, MEL} (4.21) 
ie fey fall) NG + oS} MY} (4.22) 
Ht, = -i/2 Ye {oh INR, + OER ME + oh NEN + EE MEF. 29 
Hy = ~i/ 5° {oh NE HK ME} (4.24) 


It is easy to generalize to the case of multiple shells. Every additional shell will 
have its solution expressible in the form given by Equations 4.8 and 4.23. We have a 
set of four independent equations from the boundary conditions for each interface. A 
single shell problem with two interfaces has eight independent equations from which 


2, a, a 0, BR, BRM, BON, 


In >4in >4in 1%ln > 


Similarly for a problem with N shells with N +1 interfaces, we have sets of 4N + 4 


we can obtain the sets of eight coefficients fas a 


unknown coefficients and 4N +4 equations from the boundary conditions to solve for 


them. We shall soon see that the solution for such systems are completely determined. 


4.3. Expansion Coefficients of the Scattered Wave 


We shall consider the equations of continuity of tangential E and H only. The eg and 
e, vectors will serve as the two independent tangential unit vectors. Let us consider 
the spherical interface at r = Ry. The continuity of the tangential components of the 


fields give us the following equations, valid for all points on the sphere at r = Rj: 

















E:-ey|_, = Ea: Naan (4.25) 
E; - eg = E, - eg = (4.26) 
H, - ey = Hy, - e,4 sit (4.27) 
H, - e5 Ldap, tes Hz - e9 age (4.28) 
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Similarly, the equations from the surface at r = Rg are: 








E;-e4|_,, = Es -e4|_,, +Ej-e4| (4.29) 
E,-e)|_, = Es-eo|_, +Ei-e0| (4.30) 
Hz -e4]_, = Ho-es|_, + Hi-es| (4.31) 
H,-eo|_, = Ho-es| _, +Hi-eo| (4.32) 





The expressions on both sides of the boundary condition equations above are functions 
of 6 and ¢ only. The ¢-parts are of the form e+’*. Multiplying both sides of these 
equations by e~’® and integrating with respect to ¢ over the interval ¢ € [0,27] will 
decouple all the m = 1 equations. Similarly all the m = —1 equations are decoupled. 
The expressions on both sides then reduce to being functions of # alone. Carrying out 
similar integration of both sides of these equations with respect to 6 over the interval 
6 € [0,z] could decouple the different n’s, leaving us with linear equations in sets of 
the desired unknown expansion coefficients of order n. 


Let us consider Equation 4.25: 








By-ey = So fall, MOL ep + 08H, NOY +4} sid (4.33 
n=1 

E,-eg = Do {all MEI es + 00M, NEY eg + abt, MEI, eg + OE, NEI, eg} 
=1 

(4.34) 

Now, 
Min 4] n, = —Jn (hi Ra) a Px” (cosd) ei"? (4.35) 
: a. ora P™(cos@) , 
Q)i. = Paneer ee am Vee’ ime 
And, Ni: Cp) te. IN ae (aca (alk) simp (4.36) 


Multiplying both sides by sin?6 P™(cos@)e~*”* and integrating within limits of @ € 
(0, z] and ¢ € [0, 27], and observing that 
0, (4.37) 


T ! 
and - dé sin P”” (cos@) P™(cos@) = = i < + a (4.38) 


™ d 
i d@sin?6 P™ (cos6)— P™(cos6) 
0 dé 
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the coefficients of all the a,,, terms will vanish. On simplifying both sides, we obtain 


the equation: 





a ee a ae oa 
0 Ct emacs ek _ ps eee 
+1n Gee (rj ( “) +1n kor dr (rj ( ar) = 
ld 
_ pl2yh af (1) = 
bein ( aa (rif (ur) =0 (4.39) 








Similarly, when considering Equation 4.26, only the coefficients of the a,,,, terms 
will survive. The integration procedure is identical to that used for obtaining Eqn 4.39. 


Carrying out the procedure, we arrive at: 





ayy, (Jn (1 Ri)) — ali (in(koR1)) — af (AO (ko R1)) = (4.40) 








To obtain the equations from the continuity of tangential H, we observe that by 
virtue of Equation 4.2, we can use the corresponding equations derived for the con- 
tinuity of tangential E by interchanging the a,,, and 6,,, coefficients and appending 
an extra factor of a Thus Equation 4.27 for continuity of tangential H along ey is 


obtained as: 


Gy pe ee oe ae ee 
Gin [a (as 9a 1r)) Weg Gin fla kor dr pal 2r)) ohe 
— gir [ez id (1) = 1 
A4Lin a (= re (rif (kzr)) - = 0, (4.41) 




















and Equation 4.28 gives us the continuity of tangential H along eg : 
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_ fa. oe 
ots [an (leRa) — OEM | jn (ta R1) (4.42) 
M1 b2 











Similarly the equations for the interface at R2 are obtained. Continuity of tan- 


gential E along ey gives us: 





ee ee ld 
+00 (sca rina) +H ( (“iPr ) (4.43) 


r=Rp hg a dr r=Rz2 


ld ld 
Re fez fp) — pr See aa Pr 
bein ( dr (raf (rr) > OLA (a dr (dalbor) 








And the continuity of tangential E along eg at r = Re gives us: 





(2 


ay (Jn( ko R2)) + al (AO (2 Ro)) (4.44) 


— alti (RO) (ks Re) = aap (dn( ks Re). 








The continuity of tangential H along eg at r = Ry gives: 





aj jafid,. k (yn fe2 fl a AO( hk 4A 
a+1n jig (at (Tjn(kar)) a + G4 in AD \ Bap de (r Baal 2r)) (4.45) 


_ er [E(4 4 oye nee eee, 
lS (LZ (0%) tS lea ee) 


Finally, the continuity of tangential H along eg at r = Rg gives: 
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( — HO(aa)) (4.46) 


H2 
x /€3 
Otin ( 2 jul toFa)) 


| 
Ss 
ae 
3 > 
CE 
| on 
wo 
ao 
oe 
— 
= 
wo 
ee 
iw) 
* 
Sor 
I 











Thus, for each n, and m = +1, we have a non-homogenous system of eight equa- 
tions in eight unknowns, and the system is straightforwardly solved. If we set €; = €2, 
making the shell and the core the same, then we obtain the solution to Mie scattering. 
By using the formalism discussed above, a direct solution of Mie scattering could be 
obtained by solving a simple 4 x 4 system of equations. 

In Figure 4.2 we show the result of our calculation on a glass sphere coated with 
a thin gold shell. The glass sphere (€ = 2.5) and the gold shell have a combined total 
diameter of 2 where \ = 620nm. The dielectric constant of gold at 2eV (620nm) is 
— 10.885 + 1.3482 [24]. The thickness of the shell is reduced from 0.15, to zero. When 
the gold shell is sufficiently thick, in excess of several skin depths, the optical field is 
unable to penetrate the shell. Thus the field inside the sphere will be close to zero. 
When the shell is made thinner, the optical field ‘leaks’ into the cavity formed by the 
sphere-shell system. Interestingly, the trapped photons give rise to fields inside the 
cavity that are stronger than the incident external field. The mode patterns emerge 
naturally as a result of solving the boundary value problem of the shell-core model. 
When the gold shell is made sufficiently thin, so that it can be considered a glass 
sphere only, we observe the familiar effect of lens action, the incident light getting 
focussed at the opposite end of the sphere. It is evident from these calculations that 
the screening action of a thin gold shell is remarkably efficient as far as the external 
field is concerned. The external field itself is not seriously perturbed until the shell 
thickness is made less than 0.05\. In other words one could not differentiate a solid 
sphere of gold from a glass core with a thin gold shell by observation of the external 


scattered field alone. 


4.4 Poynting Vector Calculations 


The mean intensity of energy flow at a point in an electromagnetic field with sinusoidal 


time dependence is given by the real part of what is defined as the complex Poynting 





Figure 4.2: Scattering from glass sphere coated with a thin gold shell. The 
total diameter of the core-shell scatterer is 2\. The thikness of the gold 
shell is reduced from 0.15. (top-left) to zero (bottom-right). Incident light 
is directed from the bottom of the page with polarization pointing out of 
the page. White represents larger fields. 
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vector, S: 


1 
S = 5Ex H. (4.47) 


So, the mean intensity of energy flow is given by 
a | 
s= 5Re(E <x H*). (4.48) 


Although it is only in the far field that this quantity has the simple physical picture 
of energy flow from a plane wave field, it is still instructive to study the pattern of 
energy flow in the near field. If we observe the energy flow through a spherical surface, 
concentric with our scattering object, we can study the relative energy flow from the 
wave field into the scatterer and vice versa. The knowledge of this near field energy 
flow and distribution is important for example in the understanding of radiation- 
induced interactions of adsorbed atoms/molecules on surfaces of nano-particles. 

In order to evaluate the energy flow per unit area, using the complex Poynting 
vector formalism, we need to calculate the total E and H and obtain the dot product 


of S with e,. The radial component of S is given by 
i bad bad 
5 Eo Hy — Eg He| . (4.49) 


Since the E and H fields are expressed in terms of the M and N functions, both 
of which are represented by distinct orthogonal components in the spherical polar 
coordinate sytem, it is relatively straightforward to obtain the value of the expression 
above. We show the results of such a calculation for gold clusters of radii, R = 20nm 
and R = 100nm, with 2eV (620nm) incident photons in Figure 4.3 . The observation 
point is located on circles with different radii, b, concentric with the cluster. We let b 
vary from the near field (b¥R) to the far field (b>>R). As explained earlier, by setting 
both the shell and the core equal to gold, we converge to the equivalent problem of Mie 
Scattering. It is clear that the energy flow in the near field is fairly non-intuitive. For 
certain configurations, the near field energy flow is negative, opposite to the normal 
far field scattering solution which invariably has positive energy flow outward from 
the scattering object. However, the knowledge of such electromagnetic interactions is 
vital to our understanding of certain cluster properties. Our results are in excellent 
agreement with other published works [26] on gold clusters using Mie’s theory. This 
provides us with a simple test against which to check our formalism of basis function 


expansion for the solution. 
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Figure 4.3: Poynting vector calculation for scattering from gold cluster. 
The incident light at 2.0eV is incident from 6 = 0°, with the polarization 
directed out of the paper. Scattering from two cluster sizes, R = 20nm 
and R = 100nm are shown, where b is the radius of the observation points. 
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4.5 Numerical Calculations of Total Cross-Section 


As shown in the appendix(E, we can obtain the exact expressions for the total cross- 
section and the scattering cross-section of spherically symmetric particles such as the 
single shell particle/object we have been considering. The results apply to any multi- 
shelled particles as well. The expression for the total cross section in terms of the 
scattering coefficients is : 
Oo fy? 2 

Qe = —FpRe De EE fa (aly + be RY} (4.50) 
where af, and 67, are the usual expansion coefficients of a plane wave incident along 
the z-axis with its polarization along the x-direction. The symmetry allows us to 
consider any general plane wave, but the algebra is simpler for the above choice 
of the incident plane wave. (3h and pis)h are the scattering coefficients for the 
exterior region of a single shell scatterer and should be substituted with the relevant 
coefficients for the exterior region in a given general problem. 

The net extinction or attenuation of incident energy in an electromagnetic wave 
propagating through a medium consisting of a random distribution of scatterers is 
directly proportional to the total cross section of the particles. It is assumed that the 
distribution is sufficiently dilute and random so that there is no mutual interaction 
or interference among the scattered waves from the individual particles. For many 
experimental situations, such as colloids in solution, this is a very reasonable assump- 
tion. Thus the theoretical computations on plasmon resonances in such particles can 
be verified against experimentally observed absorbance/extinction spectra. 

As another test for our solution to electromagnetic scattering from coated nanopar- 
ticles, we compare the plasmon resonance in gold colloids. Figure 4.4 shows us the 
experimentally obtained absorbance spectra of gold colloids, of approximately 6nm 
radius, in an aqueous solution. In the model calculations, we set both the shell and 
the core as gold. In the same figure we also show the calculated total cross-section of 
spherical gold particles of 6nm radius surrounded by water (n = 1.78). To do the cal- 
culations we provided the known dielectric function of gold within the range of 400nm 
to 1100nm [24]. We assumed the dielectric function of water to be a constant within 
the same spectral range. There are no adjustable parameters in this calculation. The 
only input is the particle size, which is also experimentally determined [7] by TEM 
measurements. The resonance at ~520nm is satisfactorily reproduced, including the 


aspect ratio of the peak to the satellite valley at ~450nm. The experimental graph 
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shows a peak at close to 970nm and that is attributable to water. There are small 
differences in the peak position as well as the general shape and width of the peak. 
These differences could result from several factors, including the presence of a shell 
layer surrounding the particles. Without adequate experimental verification, such 
differences could be accounted for theoretically by using our general formalism for 
the shell problem. In the following sections we shall study some of the qualitative 
aspects of the result of a dielectric shell on a core of a different dielectric. Due to 
subtle interaction between the three media, the core, the shell and the exterior, both 


the location of the peak and its shape can undergo dramatic changes. 


4.6 Test Dielectric Functions 


For calculations of plasmon resonances in single or multi-shelled particles, we need 
information on the complex dielectric functions of the various materials in question. 
Although the optical constants are accurately known for a number of materials, it 
is still not an exhaustive list. There are many materials of current research interest 
whose optical constants are not known in any great detail. In order to study model 


systems, we therefore consider the process of obtaining a ‘test dielectric function’. 
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Figure 4.4: Plasmon Resonance in 12nm Gold Colloids in Aqueous 
Solution 
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In principle, the known absorbance of a material can provide us with its complex 
dielectric function by using the Kramers-Kronig relations. In theory this demands the 
knowledge of the absorbance/extinction function for all frequencies. The conversion 
from an extinction coefficient k(v) spectrum to a refractive index spectrum is given 


by: 


An(v;) 


00) 


n(vi) — nf 
~ “Pp ih dv = (4.51) 





where P indicates that Cauchy’s principle value of the integral on the right-hand 
side is obtained. To perform a numerical Kramers-Kronig transformation, we assume 
that k(v) is given as a sequence of discrete values k;, at m frequencies v; with equal 
interval in between each v; : 

ie VE Dae as, Be as, Dag Day 

Ri, ge ies Sade, dee Sais Recs gg 

Since An(v;) has poles at each v;, we can evaluate the integral numerically by 
considering only every other point so that computations at vy = v; are avoided. Thus 
when iis an odd number (1, 3, 5,...), we perform numerical integration by considering 
v only at the even numbered frequencies and vice versa. For this method to be of any 
use, the extinction spectrum k(v) must go to zero at both ends of the entire range 
of v being considered. Such spectra are difficult to obtain experimentally. However, 
one could construct a model dielectric function by specifying a suitable extinction 
spectrum that vanishes outside the interval of interest. 

Another approach is based on obtaining an exact analytic formula for a model 
extinction function that is also specified in a closed form expression. In particular, 


we consider a double-Lorentzian function for the extinction coefficient k(v): 


kmax (y/2)° kmax (7/2) 
(v — vo)? + (7/2) (v + ¥0)? + 4/2)" 


where kmax is the maximum value of k, at frequency v9 and ¥ is the bandwidth at half 





k(v) = (4.52) 


height. It can be shown that the corresponding expression for the An(v) spectrum, 
obtained by taking the Cauchy’s principle value of the integral in Kramers-Kronig 


relations is given by: 


(vy + vo)(y/2) (vy — vo)(y/2) 
(yivo) a2)" eto) ye)? 





An(v) = n(co) + kmax . (4.53) 
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Clearly, the addition of the second term in the expression for k(v) automatically 
satisfies the conditions: k(v) = —k(—v) and n(v) = n(—v). The real and imaginary 
parts of the complex dielectric function, € = €,+7€; are related to the refractive index, 


n, and the extinction coefficient, k, by the following relations: 


Cf = a (4.54) 

€; = 2nk (4.55) 
1/2 

— {0.5 (je+d+e)} (4.56) 
1/2 

oe {0.5 (e+d-«)} (4.57) 


In Figure 4.5 we show the results of our calculations of the total cross-section of 
a gold nanoparticle of 5nm radius with a dielectric overlayer of Inm thickness. The 
dielectric function used is obtained from Equations 4.52 and 4.53 with kmaz = 1, 
y = 0.5eV and vp is made to vary between 2.2eV and 1.2eV. n(oo) is set to 1.5, a 
value typical of many organic and glassy materials that are transparent in the visible 
range of wavelengths. The important feature that emerges from this calculation is 
the appearance of a distinct second peak in the total cross-section spectrum of gold 
nanoparticles. The position and shape of this peak is distinct from that of the cor- 
responding peak in the extinction spectra (k). This illustrates the subtle interaction 
between the dielectric functions in determining the overall plasmon resonance in a 
core-shell geometry. 

In Figure 4.6 we show the calculations on a gold core of 5nm radius, with a 
dielectric shell of varying thickness. The test dielectric in this case is specified by 
kmae = 1, y = 0.5eV, n(co) = 1.5 and vg = 1.5eV. We observe significant relative 
changes in the two peaks. The plasmon resonance peak at ~ 500nm shifts towards 
longer wavelengths. Also the ratio of the peak to the valley at ~450nm increases as 
the shell thickness increases. The peak at ~750nm shifts towards longer wavelenghts 
as well. This indicates the importance of an ‘interfacial’ layer in determining the 


plasmon resonace. 


4.6.1 C60 Coated Gold Nanoparticles 


In Figure 4.7 we illustrate the subtle interaction between real dielectric media in a 
core-shell geometry. We simulate a shell of Cg on top of a gold core. The experi- 


mentally determined dielectric function for thin films of C¢q have been used for this 
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Figure 4.5: Total scattering cross section (in m”) of gold nanoparticles, 
of 5nm diameter, coated with a shell of Inm of a test dielectric whose 
refractive index, n, and the extinction, k, are also indicated. k,.: = 1, 
y = 0.5eV and 1 is varied between 2.2eV and 1.2eV. 
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Figure 4.6: Total scattering cross section (in m’) of gold nanoparticles, of 
5nm diameter, coated with a shell of varying thickness of a test dielectric 
specified by knaz = 1, y = 0.5eV and 1 = 1.5eV. The peak positions and 
shapes change as a function shell thickness. 
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Figure 4.7: Total scattering cross section of gold nanoparticles, of 5nm 
diameter, coated with a shell of Cgo. The peak at ~ 520nm shifts towards 
longer wavelengths as the thickness of the shell is increased from Inm to 
10nm. 
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calculation [32]. The striking feature in this simulation is the shift in the plasmon 
resonance peak at 520nm towards longer wavelengths as the gold core is surrounded 
by increasingly thicker shells of Ceo. The gold core is kept constant at 5nm, and the 
Ceo shell increases in thickness from lnm to 10nm. It is also interesting to note that 
the gold core makes its presence felt even when inside a relatively thick shell of C go. 
Also, the steepness of the slope and the width of the gold-plasmon peak changes as a 
function of its environment. This clearly indicates the importance of the role played 


by the medium surrounding a given nanoparticle in determining its optical response. 


4.7 Plasmon Resonance in Gold-Sulfide Systems 


There is considerable interest in the scientific community to understand the physics of 
colloids. In particular, interest in gold-coated nanoparticles has led to experimental 
observations that are not understood satisfactorily from a theoretical point of view. 
As we have seen earlier, pure gold colloid (size ~ 10nm) in an aqueous medium 
lead to a plasmon resonance peak at approximately 520nm. However, a system of 
gold-sulfide nanoparticles coated with a layer of gold exhibit remarkable dynamic 
complexity in its absorption spectra. An absorption peak, distinct from the 520nm 
resonance, undergoes characteristic non-monotonic shifts. 

In Figure 4.8 we show a typical experimental situation [7]. The formation of Au2S 
colloids is initiated by mixing together chloroauric acid (HAuCl],) and sodium sulfide 
(Nag5). This is indicated as t = 0. As the reaction proceeds, absorption spectra 
are obtained at intervals of a few minutes. Two distinct peaks are seen to evolve. 
There is the usual peak at 520nm which essentially does not undergo any shift. The 
other broad peak is at somewhat longer wavelengths and is observed to shift initially 
towards longer wavelengths but subsequently reverses its direction of shifting towards 
shorter wavelength. It has been speculated in the literature that such shifts are due to 
a combination of plasmon resonance in such colloids and resonance due to quantum 
confinement of electrons in a thin shell (of gold). It has been argued that the initial 
shift in the second absorption peak towards the red cannot be explained by a plasmon 
resonance effect. 

In Figure 4.9 we show a plasmon resonance calculation of the total cross section 
of a core-shell particle. The core has a dielectric constant of ~ 5, which is assumed 
to model gold-sulfide (Au2S) within the wavelength range of interest, and the shell is 
a thin layer of gold. The experimentally determined bulk dielectric function of gold 
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Figure 4.8: Absorbance of gold sulfide colloids coated with a shell of gold. 
The peak at ~ 520nm corresponds to gold particles. The peak at longer 
wavelengths shifts initially towards the red (1 — 4), but eventually shifts 
towards shorter wavelengths (5 — 10). 


x 10° CORE GROWTH, CONSTANT SHELL (R2-R1 = 0.5nm) 
2 





4nm 


a 
C) 


EXTINCTION(A.U.) 
a ® 


0.2 











0 
400 500 700 800 900 1000 


600 
WAVELENGTH(nm) 


Figure 4.9: Calculated Plasmon Resonance in core-shell model. The shell 
is of gold. The core is a dielectric with « = 5. The diameter of the core 
is increased from 1nm to 4nm, but the shell thickness is held constant at 
0.5nm. The peak shifts towards the red. 
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was used for this calculation [24]. The diameter of the core was increased, keeping the 
shell thickness constant. This is assumed to model the experimental situation where 
initially the formation of unstable Au2S causes the particle diameter to increase. 
However, it is argued that a ‘transition’ shell of gold is formed as a result of reduction 
of the Au-S bond on the surface by S*~ ions [15] [44]. In this ‘core-growth’ regime, the 
calculated plasmon resonance peak shifts towards the red. The peak height increases 
and this is primarily due to the larger cross sections resulting from increasing particle 
sizes. 

Figure 4.10 and Figure 4.11 are calculations for two different situations that could 
follow the initial ‘core growth’ regime. In Figure 4.10 we calculate the cross section 
when the total diameter of the core-shell composite is held constant, but the diameter 
of the core itself is made to decrease. This is the ‘diffusion mode’ which corresponds 
to the possible physical situation in which the process of Au25 formation ceases, so 
the particles cannot grow any bigger, but the process of reduction of Au2S to Au 


continues by a diffusive process through the shell, thereby increasing the net shell 
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Figure 4.10: Calculated Plasmon Resonance in core-shell model. The shell 
is gold. The core is a dielectric with « = 5. The total diameter of the core- 
shell is held constant at 20nm, but the shell increases in thickness from 
2nm to 15nm. The peak shifts towards the blue. The absorption becomes 
more gold-like. 
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Figure 4.11: Calculated Plasmon Resonance in core-shell model. The shell 
is of gold. The core is a dielectric with « = 5. The diameter of the core is 
held constant at 10nm, but the shell increases in thickness from 1.5nm to 
6nm. The peak shifts towards the blue. 


thickness. An alternative scenario is described by the calculation in Figure 4.11. In 
this ‘constant core’ mode, the core of Au2S is assumed to have stopped growing, but 
the thickness of the gold shell continues to increase as a result of the formation and 
deposition of gold on the already existing particles. In either case, the calculated 
plasmon resonance peak shifts towards the blue. It is of interest to note that the 
subtle interactions between the core, the shell and the exterior medium result in very 
non-intuitive shifts, not only in the position of the peak itself, but also on the peak 
heights and widths. Without further experimental evidence it is not possible at this 
stage to distinguish between the two possible cases, however, if one considers the 
slow growth of the Au peak during this blue shift of the Aug5, it is unlikely that the 
shell is growing sufficiently in the ‘constant-core’ model to induce such large shifts. 
The fact that one can possibly explain the characteristic blue-shifts from such exact 
calculations of plasmon resonances in the shell-model is worth further attention. 

In the above calculations, we did not make any provision for including a statistical 
distribution of particle sizes. Such distributions will lead to an overall broadening of 


the peaks. A further source of broadening can occur if the core itself is allowed to be 
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absorbing. With the availability of accurate dielectric function for Au2S as well, such 


a model can be tested theoretically against experimental observations. 
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Chapter 5 


Scattering From Non-Spherical Objects 


5.1 Introduction 


At present, solving Maxwell equations exactly for the problem of electromagnetic 
waves scattering from objects with arbitrary shapes can be tackled by finite element 
numerical methods only. The coupled partial differential equations are solved by 
finite-element methods for the vector fields E and H. Even with such powerful nu- 
merical methods there are some serious difficulties. In principle one could enclose 
the scattering object in question in a sufficiently large box with a linear or non-linear 
grid defined within it. The computed solution must satisfy the boundary conditions 
on the surface of the scattering object as well as the bounding box itself. Due to 
limitations of current computing hardware, one is restricted to a box that may not 
be made sufficiently large. This leads to an important problem in trying to obtain 
the boundary conditions on the bounding box itself. If the box is sufficiently large, 
so that the scattered solution can be assumed to be negligible at the box boundary, 
one can assume that the total field at the bounding-box boundary is equal to incident 
(plane wave) field. Since this condition is frequently difficult to satisfy, the computed 
solution will be in error. 

In this chapter, we discuss an ‘exact’ method of solution for scattering from an ob- 
ject whose boundaries do not confirm to the coordinate surfaces in a given coordinate 
system (the spherical polar coordinate system in this case). We have seen in the pre- 
vious chapter that the spherical symmetry of the shell-model allows one to decouple 
the boundary condition equations, so that one has to solve a small matrix ( 4x4 for 
Mie scattering and 8x8 for a single-shell problem) to obtain the unknown expansion 
coefficients of the scattered field. With boundaries of arbitrary shapes one has to 
solve a larger matrix, since the boundary-condition equations cannot be easily decou- 
pled. We illustrate this method of solution by considering an oblong ‘capsule’ shaped 
scatterer. Even for an object that has this relatively simple shape of a ‘capsule’, one 


previously had to rely on finite-element numerical methods for the solution. 
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5.2 The Model 


The scattering object is a ‘capsule’, composed of two hemispheres of radius R, sepa- 
rated by a cylinder of radius R and length L. The origin is chosen to coincide with the 
center of the lower hemisphere as indicated in Figure 5.1. When L— 0 the capsule 
degenerates to a sphere of radius R whose scattering solution can be obtained ex- 
actly. The axis of symmetry coincides with the z-axis. As we shall see, our choice of 
azimuthal symmetry makes the different m values decouple. The surface is therefore 
determined by three piece-wise continuous functions: The upper hemisphere section 
AB (@ € [0, 4)]), the cylindrical section BC (@ € [@,,7/2]) and the lower hemisphere 
CD (6 € [x/2,7]), where tan @, = R/L : 


r\? r L Ys 
SAB (a) @ (a) ee (a) Y et) 
Spc : R-rsné=0, (5.2) 


Sop " R -T= 0. (5.3) 


The normal direction at a given point (r,@) on the surface is obtained by taking the 
gradient of 5:n ~ VS. One of the tangential directions is chosen to be in the eg 
direction, so that the other tangent is given by n x eg. Thus the tangents and normals 


(not normalized to unit length) on the three segments can be shown to be as follows: 


n' = (r—Lcos@)e,4+ Lsin6eg 


Sap : ty; = Lsinée, —(r—Lcosé)eg (5.4) 
t, = eg 
n? = sinfe,+cos$eg 

Spo : t? = cosde, —sinbe, (5.5) 
i = e¢g 
oe 

Sep : t? = €&¢ (5.6) 
t) = eg 


The surface is approximated by choosing a set of points on the capsule, labeled 
as 1,2,3,...N. In principle, the more points we specify, the better the approximation 


to the true surface. By forcing the boundary conditions to be satisfied at the chosen 
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Figure 5.1: (a) Geometry of ‘capsule’ shaped scattering object. Two hemi- 
spheres of radius R are attached to the end of a cylinder of length L. 
Letting L— 0, the capsule degenerates to a sphere. (b) Approximating 
the azimuthally symmetrical surface by means of the discrete set of angles 
§;. r as a function of @ is specified by a piece-wise continuous function. n 
indicates the normal and t the tangent orthogonal to eg. 


points, often in a least-square sense, we can solve for the set of unknown expansion 


coefficients of the solution in the different regions. 


5.2.1 Convergence Consideration 


We have seen in Chapter 3 that the number of terms required to obtain convergence 
(to better than 1%) about the origin, for an arbitrary plane wave expansion in the 
basis functions is given by n ~ 1.4|kr|, where r is the distance from the origin and k is 
the magnitude of the complex wave-vector in that medium. This provides us a direct 
way to estimate the number of terms to include in the expansion for the scattered 
field, given the desired accuracy of the solution and the size of the scattering object. 
Clearly, the larger the scatterer, the larger the number of terms one needs to include 
to obtain the scattered solution. Such considerations are equally valid even in the 
case of the symmetrical scattering objects we encountered in the shell-model. The 
rationale follows from the following argument: Any scattered field can be expressed 


as a sum of plane waves, so if the expansion for the plane waves themselves requires a 
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certain number of terms within a region, so will the expansion of the scattered solution 
within the same region. For objects that deviate from the spherical geometry, more 


expansion terms for the scattered fields are required as well. 


5.3. Expansion in positive m only 


When representing the incident plane wave field, we used a double summation: n € 
[1,2,3,...] and m € [0,41,+2,...,+n] over the M and N functions. The ‘forced’ 
solutions of the scattered fields will have terms corresponding to each source term. 
By employing certain symmetric relations, we can express the incident field (and so 
the scattered fields as well) in a series of functions with m > 0 only. This will not only 
reduce the computation by at least 50%, but it will also avoid certain troublesome 
numerical inaccuracies with computing the m < 0 functions. The computation of 
the Associated Legendre functions for m < 0 for large |m| and n involves obtaining 
ratios of factorials like (n — m)!/(n + m)!, and these can lead to severe numerical 
inaccuracies. 

We define a set of functions Memn, Momn; Nemn and Nomn Where m > 0 only and 


the o and e subscripts refer to ‘odd’ and ‘even’ : 








Hin. <= (es) sin mo eg — pa ee) cosmé eg (5.7) 
sin 9 d 
Monn = moy( br) oe oS) cosmo eg — Lips sin mé e¢ (5.8) 
sin 
_ BRT) a. sleer dP™ (cos 8) 
Nemn = n(n+1) ras P*"(cos 9) cos mé e, + ber dp Pent) cos Mo eg 
ld P™(cos@) . 
_ mo, (rén( kr )) sin mé eg (5.9) 
n(K : 1 ue . 
Hie Se ne 12 ( ") pm (cos 6) sin mé e, + Ea tag ee) sinm@ €g 
kr kr dr 
ld P™(cos 6) 
+ mo (ren(kr)) cos m¢ eg, (5.10) 


where z,(kr) are either the spherical Bessel or the spherical Hankel functions. It can 


be verified by straightforward algebra that the following relations hold: 


Minn — MZ, = 2 Momn (5.11) 
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Moga toMy. => 2 Men (5.12) 
Nun Niuge = 2) ion (5.13) 
Nit NZ = Pile. (5.14) 
Now we have already seen that: 
m(m— my! 
Monn = (-1) Gaal! Tray M*, (5.15) 
m(ta (n—m)! 
Nege=" (oP) Gaal ae INE a (5.16) 


Since the M* and N* functions can be expressed as scalar multiples of the M and 
N functions, therefore they satisfy the diffraction equation as well. Also, since the 
m and n functions defined above are linear combinations of the M and N and their 
complex conjugates, they too are solutions of the diffraction equation. 

Consider an expansion for a p-polarized incident wave as shown in Figure 5.2. 
Both E and k are confined in the XZ plane, so that 3 = 0 (( is the azimuthal angle 
of k). Since the incident electric field, E = (F',,0, E,), do not have any component in 








P-POLARIZATION S-POLARIZATION 


Figure 5.2: Defining the relative orientation of E and k for s and p polar- 
ization geometries. The definition assumes the presence of an interface at 
the xy-plane. 
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the y-direction, we can write (since the n = 0 terms, M and N are identically zero): 


n=1m=—-n 


n=1m=0 


$+ (Be 2 yyy + Ee pny 


) Moma + (Be On + Ez bmn) N—mn} 


n=1m=0 


+Ez (digg Minn + @2 mn M-mn) + Ez (Bing Ninn + 2 mn N—mn) } 


where in the last two steps we remember to count the m = 0 term only once. Now, 


since 8 = 0 in the chosen geometry, we have 


m+1 (n 25 m)! xz 


@ ae | 
Similarly, 
! 
be L(g ier me) r 


So we have: 


fH ig i ig 
Gg. Miia Oe jc, Mi Sin 


z z 
Oo Ming dd? a My, 


—mn 


m+1 (n + m)! z 


and eee = (eal) (n = m)! mn- 
j m(ttm)!,, 
and Peas = (=1) sear mn 


Qing (Mann + (—1)"**(—1)"Mi,4) 
a®,,, (Minn — M%,,) 


mn 


Qia®,, Morn (5.18) 


an (Minn + (—1)"*4(=1)™MF,,} 
On (Minn — M*..) 


mn 


Qiaz,, Momn (5.19) 


Bian (Nmn + (—1)"(-1)" Ninn) 
Bian (Nmn + Ninn) 
OPE hss (5.20) 


12 


= bin (Nmn + Ninn) 
DUe aids (5.21) 


Substituting these relations back into the expression for the incident field, we finally 
obtain an expansion in m > 0 only. Here E, = Eo sin 6; and FE, = Eqcos6;, and the 


wave vector is k = (k,a = 7 — 6;,8 = 0). 





Be*r = 3 5. {2i (E: Ginn + Ez Gab) Momn 


(p-polarization) ae as + 2 (Ee Sinn + Ez Binn) Bemn 











Similarly, for s-polarization, with E = Eye, and k = (k,a = 7 — 9;,8 = 0), we 


obtain: 
Be™™ = Ey DS { (abn Minn + 4 mn Mmn) + (bein Nin + bin N-mn) }- 
n=1m=0 


Using the properties of sign changes in the a¥,,, and the bY,,, coefficients, we obtain 





E eikr = Eo oe >, {2 a. Memn + 22 Ce Dis} (5.23) 


(s-polarization) n=1m=0 











Now that the expansion for the s and p-polarized incident fields are known, one 
could obtain the expansion for any elliptically polarized wave from linear combinations 
of the s and p polarized wave expansions. When 2 # 0, we can always solve for the 
fields in a coordinate system for which 3 = 0, and once the field coefficients are 
determined, transform the azimuthal angle ¢ — ¢ — to obtain the solution in the 


coordinate system of interest. This amounts to simply multiplying the solutions by 
timp 
€ : 
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5.4 Boundary Condition Equations 


For the problem of scattering from a capsule, we have two regions: inside the capsule 
surface where the dielectric function is €,, and outside the surface where the dielectric 
function is €g. As in the case of scattering from spherical shells, we assume €2 to be 
a real function of the incident wave frequency. This validates the expansion of the 
incident electric field in the spherical bessel function solutions. Since the origin is 
enclosed within the surface, the solution in region 1 will have the spherical bessel 
function solutions only. 

Consider an s-polarized incident wave of unit strength ( |E,;| = 1) with the polar- 


ization oriented along the y-axis and k confined to the xz-plane (3 = 0): 


ce eee 29 
n=1m=0 
H; a a = -y 3 { rin Hae Tle at Amn a (5.25) 
2 n=1 m=0 


where Gm, = 2a¥,,, and by», = 22b¥,,, and the superscript 2 for the m and n functions 


mn? 
refer to medium 2. We assume the scattered electromagnetic fields in the two regions 


to have the following forms (expanding in the m > 0 terms only) : 


By = = a {aj/,,m Mem =e ae Tgrontt (5.26) 


n=1m=0 


Hy TT ay | > > {b}/,,m TM i ate axl mat (5.27) 
1 


n=1m=0 


E, = Sy {a2t, m2), + be nt} (5.28) 


n=1m=0 
€ 
H, = = ey . {b;,m Momn + a7 nee (5.29) 
n=1m=0 
For each point on the surface, we can write down six equations corresponding to 
the boundary conditions. As discussed before, not all of them are linearly indepen- 


dent. In fact when L = 0, corresponding to a sphere, the equations of continuity of 


the tangential fields of E and H alone will yield a linearly independent set: 


E, : ty = E, s ty + E; ‘ ty (5.30) 
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E, ‘ te => E, ‘ te + E; ‘ te (5.31) 
H, . ti = H, 7 ti + H; : ti (5.32) 
H, i. te = H, ° te + H; : to. (5.33) 


In general, we have the remaining two equations arising from the continuity of the 
normal components of D = cE and B = pH. Here n (without the subscripts) refer 


to the normal vector and not one of the field-functions : 


«, E,-n = €9E,-n+eE;-n (5.34) 


jy Hy -n = fig Hy. -n+ po H;-n. (5.35) 


Each of the terms above of the form E-t is expressible as an infinite series when 


we substitute the expression for the appropriate field expansion. Thus for example: 


E,-ti=)7> >> een Cn tite ae ti}. (5.36) 


n=1m=0 
In order to obtain these coefficients we approximate the field expansions by truncating 
the series to a finite number of terms. The number of terms to include is decided 
finally by the overall error in the total solution. The ¢-dependence in each of the 
terms of the boundary condition equations are of the form sin m@ or cos m@. On 
multiplying both sides of a boundary condition equation, say Equation 5.30, by e’”’? 
and integrating with respect to ¢ between the limits é € [0,27], only the m = m! 
terms survive. It can be easily verified that in fact both sides of these equations have 
the same ¢-dependence. For example, since t; ~ t1,e, + tigeg, and Memn +e, = 0, 
Memn + €o ~ SIN MO , Nomn + er Y Sin ME and Nomn- eg ~ sin mé@, both the left and 
right hand side expressions of Equation 5.30 will contain only sin m@ factors. Thus 
we have all the boundary condition equations above decouple for the same m. This is 
a direct consequence of the assumptions on azimuthal symmetry on the shape of the 
scattering object. For a chosen m, only n > m terms will survive due to the factors 
containing the Associated Legendre functions. So the boundary condition equations 
(after ‘m’-decoupling) will assume the form, for Equation 5.30 say: 


m+N-1 ; : ‘ , 
ye {ati, (me : t1) = aa, (m2), . t1) 7B be Era ° t1) = oe (m7 ° t,)} 
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m+N-1 


= Le ain (Minn * a) + Bron (Mian * ta) f (5.37) 


where the right hand side of the equation contains the ‘source’-terms only. Similarly, 


we can write Equation 5.31 as: 


m+N-1 
{atin (Mente) — ai (dtan  ) + tly ch eae ta) 


m+N-1 


= LE {etmn (mann ta) + bron (Minn ta) (5.38) 


Equation 5.32 becomes, cancelling the common factor —2 from the expressions for H 


in the different regions : 


m+N-1 
‘ €y 2 €2 
= fell (| a, iy) <a, | nth ts) 
n=m M1 b2 


and Equation 5.33 becomes: 


m+N-1 
: € é € 
> {el (\/ nil, t2) — ah (i Bora - tz) 
n=m M1 b2 
€ 2 


m+N-1 : 
aS {of Eig 2) + bn (yf Amt ta) (B40) 
fa b2 


Finally the boundary condition equation from the continuity of the normal component 


of D assumes the form: 


m+N-1 
pa {ati, (1 Tae ° n) — ce (€2 ies, ‘ n) 


+4, (end, -n) — O24 (egn2h-n)} 
m+N-1 


= », | Gas (egm™ mn) + bmn (e2n?!,, ° n)} ; (5.41) 


n=m 
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and the continuity of the normal component of B gives us the equation: 


ie le 1j 4 4; 2h €2 oh 
oe Ginn (fl —Nemn * n) — Gin (p2 oe : n) 


n=m fa 


1 
7 143 2h {€2_ on 
Osan (p41 TN a . n) i ae (f2 —Momn* n)| 
feat b2 


m+N~-1 ee hat 
=" fin oy Eg 8) + Bn Hay] mm) (5.42) 
n=m Ha Ha 

In these truncated boundary-condition equations, the right hand sides are com- 
pletely known, and they constitute the source terms. The left hand sides are linear ex- 
pressions in the unknowns al!,, ba, and 02", (withn €m,m+41,..,m+N-—1). 
Thus each of the equations above can be considered a linear equation in 4N unknowns, 
where N is the specified maximum number of terms to use for a specified m. The 
coefficients in these equations are functions of r and @ (the ¢ part cancels out in the 
process of decoupling the equations for the same m). So for each specified point on 
the boundary, we obtain six equations in 4N unknowns. By specifying more points on 
the boundary, we can obtain more equations. The structure of the resulting system of 
equations can be written in matrix notation which has the form shown in Figure 5.3. 

The ‘coefficient matrix’ contains blocks of 6 x 4 elements and all the blocks cor- 
responding to a given 6; are obtained from the boundary condition equations for the 
point @;. Specifying more points on the boundary or including more terms in the 
truncated expansion for the scattered fields would tend to give us increasingly better 
approximations to the exact solution but at the expense of increasing the size of the 


boundary-condition matrix. 


5.5 Solving the Boundary Condition Equations 


The procedure outlined in the previous section allows us to reduce the problem of 
electromagnetic waves scattering from an azimuthally symmetrical object, to solving 
an overdetermined set of boundary-condition equations, Ax = b. By the nature of 
the basis functions M and N (or the equivalent m and n functions), the radiation 
conditions, i.e., the boundary conditions at infinity, are automatically satisfied. This 
is of non-trivial significance since we can essentially eliminate the source of error from 
having to consider ‘approximate’ boundary conditions at the enclosing box boundary, 


which is often approximated to be at infinity. 
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Figure 5.3: Structure of the matrix derived from the boundary condition 
equations for a specified m. We obtain 6p equations from specifying p 
points on the bondary to solve for the 4N unknown expansion coefficients. 
All the blocks corresponding to say 6,, represent the six boundary con- 
dition equations for point 6,. The unknown expansion coefficients of the 
scattered fields are obtained from the over determined system of equations 
by a linear least-square method (we require 6p>6N). 
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The system is deliberately allowed to be an over determined system (number of 
rows, M > number of columns, N) since we are considering scattering objects with 
arbitrary shapes. We can solve the system in a least-square sense. It can be shown [31] 
that this essentially involves obtaining the solution of the corresponding system of 
the normal equations AA? = A™b. However, the matrix AA? will typically have 
a larger condition number (indicating greater numerical instability) as compared to 
matrix A. Instead, the preferred method of solving the linear least-square problem is 
to obtain the singular value decomposition of the boundary condition matrix, A, and 
eliminate the troublesome singular values before attempting to solve the system. 

In order to verify that this method of solution is valid, we compare the com- 
putations of the least-square solution with the exact solution in the special case of 
scattering from a sphere for which the exact solution can be obtained easily by the 
method described in the shell-problem. In Figure 5.4 we show the result of such a 
comparison. The calculations are done on a glass sphere (€ = 2.5) of radius R = 1). 
The optical field is focussed by ‘lens-action’ and the maximum field strength is ap- 
proximately 7 times the incident field strength. We have calculated the field in the 


following two cases: 


e (a) The incident wave propagates along the z-axis (6; = 0). Only m = 1 terms 


are allowed in this case. 


e (b) The angle of incidence is 45 degrees with respect to the z-axis. The incident 
wave has been expanded in a series with all allowed m. The calculation shown 


is done with m < 5 only. 


The fields are calculated using N = 12, where m+ N — 1 is the maximum value of n 
for a given m. For the least-square solution we specified 30 points on the surface of 
the sphere. For both cases (a) and (b), the solutions computed using the method of 
linear least-square are almost identical to the exact solutions as testified by the field 
calculations, validating the accuracy of this method. Although the calculation in case 
(b) involves a different set of coefficients and a larger number of terms, physically we 
expect both the solutions to be identical, with redefinition of the z-axis to coincide 
with the direction of incidence. The computations show us that indeed this is almost 
the case, the fields in case (b) are just rotated (by 45 degrees) version of the fields 
shown in case (a). The minor difference in the field patterns in the rotated case as 


compared to the non-rotated case are probably due to round-off errors involved in 
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EXACT LEAST-SQUARE 





Figure 5.4: Comparing the method of linear least-square solution of the 
boundary-condition equations with the ‘exact’ solution for scattering from 
a glass sphere (¢ = 2.5) of radius R= 1. The optical field is focussed to 
> 7 times the incident field by ‘lens-action’. The incident field is polarized 
normal to the plane of the figure. (a) 6; = 0° (top row) , (b) 6; = 45° 
(bottom row). Physically, the 45 degree incidence case should have been 
identical to the 0 degree incidence except for an overall rotation by 45 
degrees. The differences are due to round-off errors in evaluation of the 
basis functions. 
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the calculations. Since the field patterns for the ‘exact’? computation and the least- 
square computations are almost identical for both the incidence angles, the errors 
are probably not contributed from the least-square fitting procedure, but rather by 
inaccuracies (due to round-off errors) introduced by the function evaluation of the 
spherical Bessel/ Hankel functions (of complex arguments) or the Associated Legendre 
functions. 

In Figure 5.5 we show a sequence of calculations by adding the contribution from 
the next higher order term in m. The m = 0 figure is the field calculation with 
only the m = 0 term included. The m = 1 label indicates that both m = 0 and 
m = | terms are included. The field pattern becomes increasingly accurate with 
higher orders of m included in the calculations. 

In Figure 5.6 we show the contribution from each order of m separately. It is 
evident that the major contribution of the field magnitude will result from m < 4 
within the region shown in these figures. 

Having demonstrated the validity of the numerical techniques involved in these 
calculations, we now show in Figure 5.7 a sequence of calculated fields for scatterers 
that are made increasingly more oblong. The incident light in this case is made to 
approach along the z-axis. The object is elongated from a sphere by 60% (L ranges 
from 0 to 0.3A and R = 0.5). The continuity of the fields across the boundaries as 
L is increased affords visual testimony to the accuracy of the computed fields. 

In Figure 5.8 we show the calculation done on a capsule shaped scatterer with 
the incident wave approaching from an arbitrary direction. As discussed before, the 
scattering object has azimuthal symmetry about the z-axis. However the arbitrary 
incidence direction of the incident plane wave breaks the symmetry. Calculations 
are shown for a sequence of incidence angles ranging from 0° to 180°. Physically we 
expect the field patterns to vary continously as 6;, the incidence angle, changes. For 
example, when 9; is equal to 90°, the field pattern ought to be symmetrical about 
the center of the capsule, even though the origin for defining the basis functions is 
located non-symmetrically with respect to the center of the scattering object. The 


same situation applies when @ = 180°. 


5.6 The Sphere-Plane Model 


Considerable interest exists in the problem of solving for near-field scattering of elec- 


tromagnetic waves incident on a geometric structure often called the sphere-plane 
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Figure 5.5: Field calculations by including higher orders of m, from m = 0 
up to m= 5. Glass sphere (¢ = 2.5) of radius R = 1\ and 6; = 45°. Although 
m = 1 appears to be the ‘dominant’ term in this sequence, serving to 
establish the overall field distribution, the higher orders are necessary to 
obtain more accurate calculation of the field intensities. 
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m=3 m=4 m=5 


Figure 5.6: Field calculations for different orders of m, from m = 0 up to 
m = 5. Glass sphere (€ = 2.5) of radius R = 1\ and @; = 45°. The m= 1 
contribution appears to be the most ‘dominant’ term in this sequence. 
The contributions from m = 4 and m = 5 terms are going to be small. 
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Figure 5.7: Electromagnetic scattering from a ‘capsule’ scatterer. R = 0.5) 
and L is increased from zero (sphere) to 0.34. The dielectric constant of 
the scatterer is taken as 2.5 in these calculations. Light is incident from 
the bottom of the figure with the polarization directed out of the paper. 
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Figure 5.8: Electromagnetic scattering from a ‘capsule’ scatterer. The in- 
cident light approaches from different directions. From top-left to bottom- 
right the images are shown with 6; changing by 22.5° between the images. 
Although the 0° and the 180° situations are computationally different (due 
to lack of symmetry about the origin), the fields evaluate identically, which 
is expected physically. The ‘focused’ spot achieves higher intensity for 22.5° 
incidence angle. 
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model. This geometry has been used to do theoretical study of field-enhancements 
that could be responsible for the phenomena such as SERS (Surface Enhance Raman 
Spectroscopy) [3] [33] [39]. A sphere of radius R, made of a material whose dielectric 
function is €;, is separated by distance d from a plane (Figure 5.9). The half-space 
containing the sphere is composed of a material with dielectric function €2. The half- 
space not containing the sphere is composed of a dielectric with dielectric function €3. 
The incident wave propagates in region 2 and is reflected and scattered at the planar 
interface as well as the sphere boundary. Region 2 is assumed to have a non-lossy 
dielectric function, so that the plane wave expansions are strictly valid. 

The sphere plane problem has been attempted for solution by many researchers 
for nearly two decades now. The usual approach has been to obtain a solution in the 
quasistatic approximation, i.e., to solve Laplace’s equation [3] [33] [34]. However, as 
discussed before, the introduction of complex dielectric functions in the quasi-static 
approximation is artificial. Under static field configurations there cannot be any 
losses or phase delays which result from a complex dielectric function. In fact such 
equations do not satisfy Maxwell’s equations since the principle of conservation of 


energy is violated. Furthermore, such an approximation is valid only when R < 4, 








Figure 5.9: The Sphere-Plane model. The dielectric function of the sphere 
is «, and that of the plane is «;. The medium surrounding the sphere has 
non-dissipative dielectric function ¢€2. k is the propagation vector of the 
incident radiation. Such a model is often used to do theoretical studies of 
plasmon resonances on roughened surfaces. 
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even for real dielectric functions. An ‘exact’ approach for solving the sphere-plane 
model that takes exact account of retardation effects is based on solving the vector 
Helmholtz equation with the aid of a tensor Green’s function [39]. The incident 
electromagnetic wave scatters off the spherical surface as well as the planar surface. 
These scattered waves in turn scatter back and forth between the two surfaces. The 
final field is obtained by summing over the total contribution from each pass. This 
method is also computationally intensive. 

In our approach to solving the sphere-plane model we use our basis functions 
to represent the field in the three regions. Analogous to the case of the elongated- 
scatterer problem, the electromagnetic scattering problem is reduced to solving a 
system of linear equations obtained from the boundary conditions on both the surfaces. 
The boundary condition equations are easily decoupled on the sphere surface (since we 
are working in the spherical polar coordinate system) and the only task that remains 
is that of satisfying the boundary conditions on the plane boundary. However, unlike 
the problem of the elongated scatterer whose boundary can be enclosed in a finite 
volume, there is a problem in the sphere-plane problem related to the fact that the 
plane extends to infinity. The plane wave approximation implies that the fields are 
non-vanishing at infinity. So, in order to satisfy the boundary conditions on the plane, 
far away from the origin, we require an infinite number of terms in the expansion of 
the fields. 

We overcome the problem by using our knowledge of Fresnel’s formulae for reflec- 
tion/refraction of a plane wave at a planar boundary. We write the electromagnetic 
field solutions for regions 2 and 3 as the sum of two parts each. In region 2 the 
total field is that due to the incident field, E;, the reflected field (off the substrate), 
E, (assuming the sphere is absent), and the scattered field in region 2 due to the 
presence of the sphere. Similarly, in region 3 the total field solution can be thought 
of as the sum of the transmitted field, E,; (again assuming the sphere is absent) and 
the scattered field in region 3 (due to the sphere in region 2) required to satisfy the 
boundary conditions on the substrate. 

Fresnel’s formulae are obtained by satisfying the boundary conditions on a plane 
boundary for an incident plane wave [22] [12]. The problem is reformulated as follows: 
If the sphere was absent then the solution for the field in the two regions separated 
by the plane is given by Fresnel’s formulae. The incident wave gives rise to a reflected 
wave as well as a transmitted wave. In region 2 the incident and the reflected fields 


would add together to determine the total field in region 2. In region 3 only the 
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transmitted wave contributes to the field. Since all three fields are plane wave fields, 
we can express them in a series of the basis functions. However, since we already know 
the solution, such a process is not necessary. In fact we do not gain anything by doing 
this. Now we imagine bringing the sphere from infinity along the z-axis. It is going 
to be subjected to two plane-wave fields: the incident wave field and the reflected 
wave field. We know how to solve the sphere-only problem for arbitrary incidence 
angles. Therefore the two fields in region 2 can be added together to provide a single 
‘excitation’ field to which the sphere is subjected. The sphere will scatter and the 
boundary conditions on the sphere surface will determine the relation between the 
coefficients of the scattered fields in regions 1 and 2. The scattered fields in region 
2 will give rise to scattered fields in region 3 as well. Since the incident, reflected 
and refracted fields of the plane wave excitation are already made to satisfy the 
boundary conditions on the plane (Fresnel’s formulae), so now it is only the scattered 
fields in regions 2 and 3 that have to satisfy the boundary conditions on the plane. 
Since the scattered fields decay as ~ 1/1r, we need to consider points on the plane- 
boundary only far enough for the scattered fields to decrease to a small fraction of 
the incident and reflected fields added together. The remaining procedure is going 
to be identical to that followed for the elongated-scatterer problem. We set up the 
boundary conditions matrix and solve for the (six) unknown coefficients in the three 
regions. The scattered solutions in regions 2 and 3 will again be represented by 
the spherical Hankel function solutions alone, since the transmitted field due to the 
incident plane wave will be represented by the spherical Bessel function solutions and 


they will determine the asymptotic behavior of the total solution. 


5.6.1 Fresnel’s Formulae 


Consider a plane wave (polarization not specified at the moment) incident on a planar 
interface as shown in Figure 5.10. The incidence angle is 6; and the incident wavevec- 
tor is represented as k;. 6, and 6; correspond to the reflected and refracted angles 
and k, and k, represent the reflected and transmitted wavevectors. If the wavevector 
is complex in a passive medium then of course the wave will be attenuated. Thus 
the incident, reflected and refracted (plane wave) fields in the two media takes the 


following forms: 





Ep= Ei et(hir-) H; = —_V xE; (5.43) 


tw p2 
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Figure 5.10: Plane wave incident at a planar interface. Reflection and 
refraction occurs. The plane is located at z = d. The phase must satisfy 
spatial invariance on the plane. 








E, = Ej ellkrr-wetsr) H, = Vv xE, (5.44) 
E, => He ei(Ker—w tte) H, => aa x E, (5.45) 


The fields have to satisfy definite boundary conditions at the interface. Spatial in- 
variance on the plane (by symmetry, the same relation must hold for every point on 
the plane) requires that at all times the phase factors k;-r, k,-r+6, and k,-r+ 6; are 
equal. Otherwise the field relations would depend on the position on the plane, which 
is physically disallowed. Since all points are equally valid on the plane, we consider 
the point of intersection of the z-axis with the planar boundary. We can immediately 


write (for this chosen point): 


k;-r = k;-(-de,) = k,d cos 6; (5.46) 
k,-r = k,-(-—de,) =—k,d cos@; (5.47) 
k; oD. k, (-—d e,) ml ky d cos 6. (5.48) 


Using the equality of the phases, we immediately arrive at the expression for 4, : 


6, = 2k; d cos; (5.49) 
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Using Snell’s Law: 
ij 


— sin 6, Or (5.50) 
ky 


cos; = 4/1 —sin? 6. (5.51) 


If &;/k, is a complex number then we obtain a complex number for cos 6; from the 


sin 6; 


equation above. This in turn will yield a complex number for the phase in region 3. 


So, we can write the expression for 6; as: 
6, = k; d cos; — ky d cos &. (5.52) 


Having determined the phases, we can compute the field in the two regions (2 and 
3) that satisfy the boundary conditions. The relation between the amplitudes Ej, E% 
and Ej is given by the Fresnel’s relations [22]. For s-polarization (E directed out of 
the paper in Figure 5.10) with n = \/ep: 

















Ei = 2n cos 8; (5.53) 
FG nz cos 9; + u2.4/n3 — n3 sin? 6; , 

. He Dg 2 ee aN. 
Et nq cos 6; u2./n3 ng sin’ 6; (5.54) 
Eb nz cos 9; + H2.4/n3 — n2 sin? 0, 


The corresponding expressions for p-polarization (H directed out of the plane of the 























figure) are : 
Ei 2ng nz cos 0; 
Ei” tan 2 2 ain? (5.55) 
0 ia 13 COS 6; + nav/nd — nj sin’ 6; 
2 gy 2 ; 2 2 ats 2): 
E _ an3 COs 6; — navn — nz sin* 6; — 
Fi ua? 6 2 tas. (38) 
0 fa t3 COS i+ nay/nd — 3 sin* 6; 


With the knowledge of the plane wave fields satisfying the boundary conditions on 
the substrate, we can readily handle the case of satisfying the boundary conditions 
for the scattered fields due to the sphere. 


5.6.2 The Boundary Condition Matrix 


In Figure 5.11 we show the schematic of the matrix derived from the boundary con- 


dition equations. As a variation in the numerical techniques for solving this problem, 
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we shall try to obtain the solution by solving a square matrix derived from the bound- 
ary condition equations. This is suitable when the lengths involved (R and d) are 
much smaller than the wavelength of the incident radiation. When these lengths be- 
come comparable to the wavelength, higher precision floating point calculations are 
necessary for preventing excessive round-off errors. We try to obtain the solution 
in a truncated expansion with the maximum value of n = m+ N— 1. Thus the 
matrix we are required to solve will have dimension 6.N x 6.NV. Since we require that 
6N =4N +6Nous ( a square matrix ), where N,,, is the number of specified points 
on the substrate, N must be divisible by 3. 

The first 4.N equations correspond to the boundary conditions on the sphere sur- 
face. As in the case of a simple isolated spherical scatterer, these equations decouple 
with respect to different n’s and so the first 4.N rows consists of ‘diagonal’ blocks 
of 4 x 6. The remaining elements in this upper section of the matrix are zero. In 
other words, we are combining all the individual 4 x 4 equations for a sphere and 
then attempting to solve them all at once. Within these blocks only the entries cor- 
responding to the coefficients in regions 1 and 2 will be non-zero. We require 6.V.u» 
equations for the substrate corresponding to the N,. specified points. For the sub- 
strate equations, the elements corresponding to the coefficients of region 1 will be 
zero. Convergence of the solution can be checked by either increasing the number 
of points on the substrate or changing the location of the points themselves. Either 
way, the solution may not be perturbed by more than the error limits defined for the 


problem at hand. 


5.6.3 Field on the z-axis 


On the z-axis, i.e, when 6 = 0 or 7, evaluation of the total electric field is made simpler 
because many of the terms in the expression for the total electric field vanishes. When 
k lies entirely in the xz-plane, for s-polarization (E directed along the y-axis), E can 


be written as a linear combination of Men, and Ngmn functions. Now, 


dP™ (cos 9) 


P™ (cos 6) 
dé 


Memn = —mMz,(kr) sin md e€9 — zn(kr) cosmdé eg. (5.57) 


Clearly, the first term is going to be zero on the z-axis, for any m #4 +1. For m > 0 
the P™(a = +1) = 0, and for m = 0, the factor m ensures that the first term is 
always zero on the z-axis. As we have seen before, on the z-axis, the second term 


containing the derivative of the Associated Legendre function wll be zero for m # 1. 
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Figure 5.11: Structure of the matrix derived for the sphere-plane model 
from the boundary condition equations for a specified m. We obtain 6N 
equations in 6N unknowns by specifying N/3 points on the substrate. The 
equations for the sphere boundary (indicated by the dotted box) are exact. 
The matrix above is shown for n = 6 and so the number of points specified 
on the substrate, N,,.,, = 2. 
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Similarly, by examining 





k 1 mn 0 
Nomn = 2(n+ i) P"(cos 8) sin m@ e, + Te (rag( br) es) sin M@ eg 
se P™ (cos 9) 
+ mo, (ren( kr) cosm@ eg, (5.58) 


we notice that the first term is going to be zero for any m because of the product 
P™(+1) sinm@. The second and third terms are going to be non-zero on the z-axis 
only for m = 1. Thus, for the s-polarization case the total electric field on the axis 
will not have any radial component. This is physically reasonable since the incident 
electric field does not have any radial component on the z-axis either. The total 
electric field on the axis can therefore be obtained by solving the m = 1 matrix only. 

In the case of p-polarization, the scattered electric fields are expressed in terms of 


the Momn and Nem, functions. Now, 





o P™(cos 8) dP™(cos@) . 
Hives. <— min(kr)— cosm@ eg — 2n(kr)—— sinm@ eg. (5.59) 
and, 
Nemn = (nt 1) EY) pm (cos 6) cosmé e, + Te (reg( br) es) cos md eg 
ld P™(cos@) . 
— mo (rén(kr)) sin m@ eg, (5.60) 


As for the Mym, functions, only the m = 1 terms are non-zero on the z-axis. For the 
Nemn functions, the first term is going to be non-zero for 6 = 0,7 for m = 0 only. 
The eg and eg components are going to be non-zero only for the m = 1 case. Thus, 
in the case of p-polarization, the total radial electric field will be contributed by only 
the m = 0 solution, and the total electric field in the direction orthogonal to e, will 
be contributed by the m = 1 solution only. 

This turns out to be a great simplification in the computation of plasmon reso- 
nances in a sphere-plane geometry. Since the electric field at various points in the 
same medium are linearly related to each other, evaluation of the field on the axis to 
determine plasmon resonances in a given geometrical configuration is a sound choice 


from the point of view of numerical efficiency. 


5.6.4 Plasmon Resonances in The Sphere-Plane Model 


We calculate the plasmon resonances in the sphere-plane model for three systems: (a) 


Silver sphere on silver substrate , (b) Gold sphere on gold substrate and (c) Copper 


93 


sphere on copper substrate. We use the experimentally determined bulk dielectric 
functions for gold, silver and copper [24]. The incident light is s-polarized and the 
incidence angle is 45°. Calculations are done with a sphere of radius R = 30nm, the 
sphere-plane separation d = 5nm and the ‘observation’ location is at the intersection 
point of the substrate and the z-axis. We have used n = 12, and tests with higher 
number of basis functions do not alter the results, confirming that the solution had 
converged. We choose points on the surface at fixed angular intervals. The boundary 
conditions are made to satisfy exactly at these points. If convergence is reached then 
the final solution should not depend on the particular choice of these chosen boundary 
points. 

We compare the results of our computations with the experimental data of Berndt 
et.al. [10]. The experimental data shows the spectrum of light emission from a tip- 
substrate geometry in the ‘field-emission’ regime (the tip-substrate voltage difference 
is approximately > 100Volts. Presumably, the field-emitted electrons would excite the 
plasmon resonances in the tip-substrate geometry. The theoretically obtained data 
on the plasmon resonances is weighted by the detector sensitivity graph of Berndt 


et.al. [10]. 
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Figure 5.12: Plasmon resonance in sphere-plane model. A Gold sphere 
of radius 30nm located 5nm from a gold substrate. the sphere are gold. 
(a) Theoretical calculation of plasmon resonance based on experimentally 
obtained bulk dielectric function for gold; the broken line shows the as- 
sumed detector response in the photoemission experiment in (c). (b) 
Multiplying the calculated spectrum in (a) by the detector response. (c) 
Experimentally obtained spectrum of photon emission from a tip-substrate 
geometry in the ‘field-emission’ regime. 
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Figure 5.13: Plasmon resonance in sphere-plane model. A Copper sphere 
of radius 30nm located 5nm from a copper substrate. the sphere are 
copper. (a) Theoretical calculation of plasmon resonance based on ex- 
perimentally obtained bulk dielectric function for copper; the broken line 
shows the assumed detector response in the photoemission experiment in 
(c). (b) Multiplying the calculated spectrum in (a) by the detector re- 
sponse. (c) Experimentally obtained spectrum of photon emission from a 
tip-substrate geometry in the ‘field-emission’ regime. 


The peaks in the computed plasmon resonance spectra correspond very closely 
to the peaks in the photon emission spectra in the field-emission experiments. The 
differences observed could be due to several reasons, including contributions from 
other processes, inaccuracy in the assumed detector response function, artifacts of 
the detection system (such as observing higher order diffraction peaks) and round-off 
errors in the computation itself. However, the strong correlation suggests that the 
light emission process in the ‘field-emission’ regime is due to a plasmon-resonance 
effect. 

The plasmon-resonance peak in the case of a silver sphere over a silver substrate 
exhibit large field enhancements. In the case of this calculation with s-polarized 
excitation, the enhancement is by a factor of ~ 10. This implies that the effective 
intensity of the optical field is enhanced by a factor of ~ 100. This is in agreement 
with the quasi-static calculations done by Arvind and Metiu [3]. Processes such as 
surface enhanced Raman scattering of molecules that depend on the square of the field 
intensity can therefore exhibit considerable enhancements (> 10+) due to nanometric 


irregularities on substrate surfaces. 
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Figure 5.14: Plasmon resonance in sphere-plane model. A silver sphere 
of radius 30nm located 5nm from a silver substrate. The sphere is silver. 
(a) Theoretical calculation of plasmon resonance based on experimentally 
obtained bulk dielectric function for silver; the broken line shows the as- 
sumed detector response in the photoemission experiment in (c). (b) 
Multiplying the calculated spectrum in (a) by the detector response. (c) 
Experimentally obtained spectrum of photon emission from a tip-substrate 
geometry in the ‘field-emission’ regime. Notice the field enhancement at 
resonance (by a factor of ~ 10). 
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Figure 5.15: Comparison of field enhancements with a silver sphere on a 
copper substrate and a copper sphere on a silver substrate. In this case, 
it is clear that the enhancements are contributed more by the nature of 
the sphere as compared to the substrate. 
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The large field enhancements in the case of a silver sphere over a silver substrate 
leads to the following question: Which of the two, the sphere or the plane, is more 
important in this enhacement process. It has been customary to assume in the existing 
literature to attribute such field enhancements to the substrate. It has been argued 
that the role of the particle (the sphere or a tip) has been merely to break the 
symmetry (translational invariance) of the plane, so that the surface plasmons can 
couple to the radiation field by allowing conservation of momentum to take place. To 
answer such a question we compare the results of our calculation in a sphere-plane 
geometry in (a) a copper sphere on a silver substrate, and (b) a silver sphere on a 
copper substrate. The results are shown in Figure 5.15. The incident radiation is 
again s-polarized and the incidence angle is 45°. The sphere diameter is 30nm and 
the sphere-plane separation is 5nm and the field observation point is at the point 
of intersection of the z-axis and the substrate. It is clear from this calculation that 
in this case at least, it is the partilcle (the plasmon resonance of the particle) that 
is more important than the substrate. However it is a truly coupled phenomena. 
For example, the enhancement in the case of a silver sphere on a copper substrate, 
although large (~ 6), is still smaller than the the enhancements with a silver sphere 


on a silver substrate. 
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Chapter 6 


Conclusion 


In this thesis we have been developing the technique of vector basis function solution 
of Maxwell’s equations. In all the different ‘example’ cases in which this approach 
has been applied, such as electromagnetic scattering from concentric shells or the 
calculation of plasmon resonance in a sphere-plane geometry, the incident electro- 
magnetic field has been assumed to be plane waves. The increasing sophistication 
of experimental techniques has brought the domain of ultrashort (sub-picosecond) 
optical excitations into many different contemporary experiments. The plane wave 
approximation may not be valid in these situations for several reasons. The dielec- 
tric relaxation times (which indirectly determine the dielectric functions) are often 
much longer than the duration of the optical excitation itself. We are speaking of a 
truly transient response in the sub-picosecond time scales. Furthermore, the extent 
of the optical ‘wave-train’ may often be shorter than the dimension of the scattering 
objects themselves. Also, we have not considered the problem of solving the Maxwell 
equations in the presence of sources, namely charges and current densities. Although 
such general problems are certainly not of mere academic interest, the mathematical 
and computational tools available to solve these problems effectively are still far from 


being satisfactory at the present. 


6.1 Transient Excitation 


The basic problem of transient response can in principle be handled by the method 
of Fourier decomposition of the source terms. Typically an experimental situation 
consists of a stream of ultrashort pulses, whose Fourier decomposition is determinable. 
One could in principle solve the problem for each of these component excitations up 
to some limiting order and obtain the final sum. In principle this can be readily 
accomplished for a linear medium. In nonlinear media the summation process will 


not be straightforward. 
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It is the author’s belief that a more elegant approach of solving the transient 
problem is possible by further research to find a suitable ‘Laplace-like’ transform 
(Generalized Fourier Transform) to represent the pulsed excitation. Since the ‘basis’ 
functions form a complete set, therefore the basis set satisfies the property of ‘con- 
vergence in the mean’ (this is in fact a weaker property than completeness). This 
however validates the process of interchanging the order of an infinite summation 
and an integration to obtain the generalized Fourier/Laplace transforms [21]. So the 
transform of an infinite sum is equal to the infinite sum of the transforms of each 
term. The method should be similar to that of solving for transients in electrical 
networks consisting of linear, lumped and time-invariant elements. 

The importance of the transient problem lies in its potential to elucidate ex- 
perimentally, the dynamical properties of nanometric systems (such as the scaling 
of dielectric functions with particle sizes in the nanometer regime, for example) by 
comparison with theoretical calculations. With our increasing dependence on the 
scanning probe techniques to understand basic physical phenomena on surfaces, the 


importance of the transient problem is clearly established. 


6.2 Related Problems 


In this thesis we have developed the method of basis function solution in the spherical 
polar coordinate system. We have demonstrated the feasibility of using this basis set 
for solving boundary value problems in which the boundaries may not conform to the 
coordinate surfaces of the spherical polar coordinate system (the problem of scattering 
from the ‘capsule’-shaped scatterer or the sphere-plane model for example). 

In principle, we could solve the problem of plasmon resonances in what is called the 
tip-substrate geometry using the same techniques. Such a geometrical configuration 
is of importance in the modeling and understanding of electromagnetic interactions 
at the tunneling junction of a scanning tunneling microscope. Such an understanding 
is important in many of the other scanning probe microscopies as well, such as AFM 
or SNOM etc., when light is coupled into the tip-substrate junction. 

If we could allow the sphere in the sphere-plane model to assume the shape of 
an elongated scatterer that is sufficiently long, in principle we could model the tip- 
substrate problem. However, such an approach might require more sophisticated 
numerical techniques to solve the boundary condition equations than those employed 


in this thesis. The elements in the boundary condition matrix contains products of 
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spherical Bessel or spherical Hankel functions of the first kind and the Associated 
Legendre functions. The columns in these matrices contain these functions with the 
index n increasing from left to right. Now these functions scale as a function of n and 


z as follows: 
Ald 


DEN SS ee eine) (6.1) 





ee 129-5 2(Ie= 1) 6.2) 


gntl 





hn(z) = Jn(z) +2nn(z). (6.3) 


This kind of dependence makes the matrix highly ill-conditioned. The problem is 
similar (if not worse) to that which is encountered when trying to invert matrices 
of the type known as Vandermonde matrices [31]. These arise in connection with 
problems of polynomial fitting. When we need to include more terms (to achieve a 
better approximation), the matrix could become unduly lopsided from left to right. 
This will invariably lead to numerical instabilities in the solution when the ratio 
of the largest to the smallest elements exceeds the available machine precision. For 
scattering geometries that deviate significantly from the spherical coordinate surfaces, 
one would need to include more terms. This will lead to a more ill-conditioned matrix. 

One way to overcome such a problem may be to use higher precision floating point 
computations. Software libraries exist in Netlib and elsewhere to allow one to perform 
quadruple precision floating point calculations or even arbitrary precision calculations. 
Trade-off will have to be made with respect to the desired accuracy of the solution 
and the increased computation loads resulting from the higher precision data types. 
Probably it is only necessary to perform the process of solving the matrix in higher 
precision, since it is the most critical step in determining the overall round-off errors. 

Another alternative could be to use a coordinate system whose coordinate surfaces 
are closer to the boundaries of the problem in question. The ‘recipe’ for obtaining the 
basis set without any reference to a particular coordinate system has been outlined in 
Chapter 2. That such a set will be complete has also been proved without reference to 
a particular coordinate system. For example, in the tip-substrate geometry one could 
use the prolate spheroidal coordinate system. The main handicap at present is the 
lack of adequate analytical treatments of the eigenfunctions of the scalar Helmholtz 


equation in the prolate spheroidal coordinate system. This has been the primary 
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justification in this thesis to develop the subject in the spherical polar coordinate 
system. The eigenfunctions in the spherical polar coordinate system are composed of 
functions like the spherical Bessel and Hankel functions and the Associated Legendre 
functions. The supporting literature for these functions is very satisfactory. However, 
there are a few isolated works that have a fairly comprehensive account on the subject 
of spheroidal functions and these will provide the necessary groundwork on which 
to build the subject of basis function expansion in the prolate or oblate spheroidal 
coordinate systems [16] [6]. 

Other related problems that are of interest are the solution of electromagnetic 
scattering from roughened surfaces (the substrate has a hemispherical bump), the 
problem of electromagnetic interaction among clusters in nanometer proximity to 
each other (many coherent scatterers). There is also the class of problems related to 
magnetic shielding. The present apparatus of vector basis functions can allow exact 


solutions to a shielding problems, be it electric or magnetic shielding [27]. 


6.3. Antenna Problem 


We now address the important problem of solving the Maxwell equations in the 
presence of sources. In the process of developing the general formalism for scattering- 
like problems (homogenous problem), we created a complete basis consisting of the L, 
M and N functions. However, with zero-divergence plane waves and in the absence 
of sources, the L functions did not seem to play any role other than its contribution 
to make the total set attain the property of mathematical completeness. All our 
scattering solutions were based on the M and N functions only. As it turns out it is 
the L functions that provide the missing link to solving the inhomogenous Maxwell’s 
equations. 

Maxwell equations in a linear, homogenous, isotropic and time-invariant medium 


(that is € and yw are not functions of r or t), are of the form: 


V-E(r,t) = eet) (6.4) 
V-H(r,t) = 0 (6.5) 
Vx E(r,t) = goed) (6.6) 


ot 
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Vx H(r,t) = J(r,t)+e (6.7) 


Taking divergence of both sides of Equation 6.7, since divergence of curl is zero, we 


obtain the continuity equation for electric charges: 





0 
V-I+en(V-E) = V-(VxH)=0 (6.8) 
Op 
V-J+— = é 

Or, - At 0 (6.9) 

Taking curl of both sides of Equation 6.6 and Equation 6.7, we obtain: 

OJ CE 
2 
VxVxH = Vi her (6.11) 
Or, assuming sinusoidal time dependence e~*", so that + = —iw, and ae = —w’, we 
obtain: 

(twuJ) -VxVxE+h’E = 0 (6.12) 
(VxJ)-VxVxH+k°H = 0 (6.13) 


where k? = w?1e. Now, Equation 6.12 would resemble the form of the diffraction 
equation (Eqn 2.18), if 


V(V -E) = Me = wydJ 
€ 


Or; Vp = wed 
Or, V’p = twpueV-JI 
Or, V7p = iwpe(iwp) (From Eqn 6.9) 
Or, Vptkh’p = 0. (6.14) 


Similarly, Equation 6.13 would resemble the form of the diffraction equation (Eqn 2.18), 
if 


Veh V(V +H) =o (Since V -H = 0). (6.15) 


From Vp = wyeJ, obtaining curl of both sides of the equation and noting that the 


curl of gradient is zero, we have V x J = 0. So, if p(r,w) satisfies the scalar Helmholtz 
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equation, then both E and H will satisfy the diffraction equation. In other words, 
both E and H can be expressed in a series of the {L,,M,,N,} functions and p can 


be expressed in a series of the scalar solutions {w,,}. Thus: 


p(rj2) = dave (6.16) 
1 





J(r,w) = a oe Vin = on L, (6.17) 
E(r,w) = So {ci Ln +a5M, + 05N,} (6.18) 
H(r,w) = D> {cL +ahM, +02N,}, (6.19) 


where the coefficients of the expansion can be determined by matching the boundary 
conditions. This is the solution to the famous antenna problem. We have the 
additional boundary condition stating the convervation of charge, Normal component 


of J is continous at a boundary : 


Thus, given a certain antenna geometry, the charge density can be obtained by solving 
the scalar boundary value problem for p. This in turn will give the solution for 
the current densities on the boundary of the antenna. Alternatively, the boundary 
conditions on J can be invoked. The tangential component of J enter the boundary 
conditions on E and H. The E and H fields can therefore be obtained by matching 
the boundary conditions between the assumed solutions across the boundary. As it 
turns out, p = 0 does indeed satisfy the conditions of this problem, and this subset 
is what we were concerned with, for solution to the scattering problems. Thus our 
original nomenclature of calling Equation 2.18 as the ‘diffraction equation’ is actually 
a misnomer as the class of electromagnetic problems it represents is more general 
than that required to study diffraction and scattering only. 

The vector basis function solution to the inhomogenous Maxwell’s equations will 
open up the path to handling complicated radiation problems and their near field 
behavior. Such problems are extremely difficult to handle with existing techniques 
to solve the problem. To the best of our knowledge, this is the first time the antenna 
problem has been shown to be solvable as a boundary value problem using a complete 


set of basis functions. 
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6.4 Final Remarks 


We have developed the technique of basis function expansion in the spherical polar 
coordinate system and solved a few representative problems. The problem of plasmon 
resonance in colloids, diffraction from an elongated scatterer and the sphere-plane 
model have been solved using one and the same technique. Even the inhomogenous 
Maxwell equations are solved (in principle) using the same techniques. Although the 
method of expanding solutions in ingeneously created functions has been around for 
some time, this thesis establishes the aspect of applicability of this method to very 


general problems. 
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Appendix A 


Helmholtz Equation in Spherical Polar 
Coordinates 


In spherical polar coordinate system, the scalar wave equation is readily solved by 


the method of separation of variables. Thus, given the wave equation 
V7u(r) + k’u(r) = 0, (A.1) 


we assume the existence of solutions of the form u(r, 6, ¢) = R(r) P(@) ®(¢). In 


spherical polar coordinates, the Laplacian operator has the form: 


tate NOY toe Oe OO OE AD 
Veana|> pene ag oo) oe kane oes | | 


On substitution of u = RP® into Equation A.1, and dividing by RP®, we obtain: 


1 1 Of. AP ®) 1 @(P®) ou 
ma = POr? sin @ E (sine oo ss sin@ 0¢? He ee a) 














A.1 Radial Equation 


Introducing the separation constant A 





1 1 Of. ,(P®) 1 (Po) 
a f) —— A.4 
P®siné E (sn Oo a sind 0¢? ae =) 
we obtain the separated radial equation: 
ar aa, 


If k? £0, we let r = p/k and the equation becomes 


PR 2dR x 
A ee), A. 
raarradl >) ° oe 
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If we let R = S/,/p, then the radial equation can finally be written in the form: 


OS $Vds BA 

— + —-— 1-—|S= A. 

tea 5) " a 
where B = ,/A+ . This is Bessel’s Equation and the solutions are the well known 


Bessel and Hankel functions S(p) = AJa(p) + BHY)(p). Equivalently, 
1 1 
R(r) = A—=Jp(kr) + B= HO?) (kr A8 

where A and B are arbitrary constants and the 1 and 2 in a? indicate the Hankel 
functions of the first and second kind respectively. When k? = 0, corresponding to 
dc or quasistatic conditions, we obtain Laplace’s equation and the radial equation 
becomes much simpler: 

1 al? A 

See —~—R=0. A. 

= gpa bt) att 0 (A.9) 
Assuming solution of the form R = r*, we obtain a quadratic equation for the roots 
of a, so that ag = +(—1 +V/1+ 4X). So the general solution to the radial equation 
for k = 0 will be of the form: 


Rir) = Ar®+ + Bre. (A.10) 


A.2 Angular Equations 


From Equation A.4 we immediately obtain: 





1 od dP 1 1€® 
—— — | sin 8 = 0. A.11 
P sind dO (ss a) + manasa ” oe 
The ¢-equation separates if we let 
1d 2 
so that 
0(¢) =Ce'™? + De"? (A.13) 


If the full range of ¢ € [0,27] is allowed in a given problem, then the requirement of 
single-valued solutions imply that m can only assume integer values me {0,+1,+2,+3,...}. 
The corresponding 6-equation becomes 


@P cos@dP m? 
= — 0). A.14 
d6? Pane atl je ( ) 
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On changing variables to z = cos 6, the 6-equation transforms immediately to Legendre’s 


equation: 





dz? 1—2x2 dx 1-2 ~ [~g? 


@P 2% dP 1 2 
ae (. a ) Po, (A.15) 


with solutions of the form P(x) = EP™(x) + FQ™(x), where \ = n(n +1). The 
solutions are the well known Legendre (and Associated Legendre) functions P(x). 
The second independent solution to Legendre’s equation are the Q(x) which are 
singular on the z-axis. A problem in which the z-axis is not included must have these 


solutions included in the expression for the general solution. 


A.3 Ejigenfunctions of the Helmholtz Equation 


We substitute \ = n(n +1) in the expression for the solution to the radial equations. 
Since 6? = A+ iy we obtain 8 = n+ 3, and the solution to the radial equations 
assume the form: 


1 i 


AG) = ARP Ae) PSOE 
(r) TE ng (hr) + ae eat r) 
= A'jn(kr) + BRO) (kr) (A.16) 
where 
3 via 
In(p) = Dp nt $l) (A.17) 


(12)(,) — /2 p72) 
are the spherical Hankel functions of the first and second kind. Therefore, in a region 
that includes the origin and so the z-axis as well, and where ¢ can range over the full 
azimuthal range [0,27], the eigenfunctions to the scalar Helmholtz equation assume 


the following form: 
u(r, 0,6) ~ jn(kr)P™ (cos O)e"?, (A.19) 


with n € {0,1,2,...} and m € {0,+1,+2,...,+n}. For a region that excludes the 
origin, but not the z-axis, the spherical Hankel functions are to be used instead of 


the spherical Bessel functions. 
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Appendix B 


Spherical Bessel Functions 


It has been shown that the radial part of the solution to the scalar Helmholtz equation 


consists of the spherical Bessel and spherical Hankel functions: 


B.1) 


=> 
Ss we 
eos 
aN — 
8 8 
—— —— 
| | 
ao 
+ 
Nik 
— 
8 
— 


Chet nye) 


(B.2) 
Te Inga (a) + fens (2): (B.3) 


Here n,,(x) are the spherical Neumann functions, derived from the Neuman functions, 


Ni({x), as indicated above. The argument z can be complex in general. In the defi- 
nition of the spherical Hankel functions, the first and second kinds correspond to the 


positive and negative signs respectively. Explicitly, 


fla) = aren SS LVF) on 


d B.4 
A mi(2n + 2m + 1)! i (B.4) 





1 Oa Di ge 
fina) = = Srgeel be x ‘ (B.5) 


= mi(n —m)! 





Clearly, when x > 0, jn(z) ~ x” — 0 for n £ 0. So only jo(x) is non-zero as x — 0. 
The spherical Neumann functions, n,(z) ~ 1/x"t! — oo as x > 0. The convergence 
is slow for large arguments, and it is computationally desirable to expand in a series 
of inverse powers of x. This can in fact be done and the resulting series is actually 


expressible in a finite number of terms [42]. Thus we have “ > 0): 














ee ee (FAD (7, + byt n —(n-k+1)(n 4 )! 
Gal®) = ae € ra k(n = k)\(Qa)* sods a —k)! (2 x) , (B.6) 
We) = aig Ee _ an Ba 


WED os oe n+k)! 
a lea ios Hien ae 
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Thus, unlike the normal Bessel and Hankel functions, the spherical counterparts are 
amenable for exact numerical computations that does not involve truncation errors 
resulting from limiting the number of terms used in the computation of the series 


expansion for these functions. As  — oo, we have the following asymptotic behavior: 








A(x) ~ =(-aytte (B9) 
x 

(2) i ‘n+l —ix 

AS (a) ~ -1""e (B.10) 
x 
1 1 

gala) es — cosa — "= 1) (B.11) 

na(t) ~ ~sin(e@—"**n). (B.12) 
x 


A function f(ka — wt) represents a wave travelling in the direction of increasing 
x. In other words, as t increases, the constant-phase points correspond to increas- 
ing x values. Similarly f(kx + wt) represents a wave travelling in the direction of 
decreasing x. Now with the convention of allowing the time dependence of the ex- 
citation to be of the form e~"* (as opposed to e'*), a wave ~ eet = et(kr—wt) 
represents an outgoing (increasing r) spherical wave. Thus with the convention of 
time-dependence ~ e~“*, we immediately recognize that the functions that describe 
the correct asymptotic behavior are the spherical Hankel functions of the first kind. 
The spherical Hankel functions of the second kind represent incoming waves and are 
not appropriate for expressing solutions to scattering problems. The spherical Bessel 
and spherical Hankel functions of the first kind are two independent sets of functions 


in which our scattering solutions can be expressed. 


B.1 Recurrence Relations 


The recurrence relations of the spherical Bessel and spherical Hankel functions follow 
directly from the corresponding recurrence relations of the normal Bessel and Hankel 
functions. If Z,,(a) represent either the Bessel function or the Hankel functions, and 


z,(x) represent the corresponding spherical counterparts, then, we have: 


Zy-s(2) + Zpsr(2) = P2(2). (B.13) 


Multiplying both sides by \/m/(2z) and letting p = n + 4, we obtain : 





7 7 2n+1 [x 
Fie oe ee = 4/—Z,11 : 





2n+1 
Pade) + Znti(£) = an ee) 














Similarly, by using the recurrence relation: 








als) Vip, (2) — Zyl 
we obtain 
wale) - = 1 [n 2n—1(a) — (rm + 1)2n4a(2)| 











In addition we have the following recurrence relations: 

















d 41 n+1 
ae [2 zn(2)| = 2" 2n-1(2) 
ine ee =n 
ae zn(2)| = —2"Zn41(z) 
Now, 
1 (req(kr)) _ <sa( kr) 4 oll) 
= NS negli) = (nt Henealbr) 
on re (Zn—1(kr) + 2nti(kr)). 
Thus, 
ld 
wp (Peal) = Sy ((n + Len-a(hr) — menea(hr)) 
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(B.14) 


(B.15) 


(B.16) 


(B.17) 


(B.18) 


(B.19) 


(B.20) 


(B.21) 


(B.22) 
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B.2 Special Integrals of Spherical Bessel Functions 


A fundamental integral involved in obtaining the expansion coefficients of an incident 
plane wave in the L, M, N basis, is the following: 


- J,(at) J, (at) “ sea) (B.23) 


- vy? = 2 





Substituting the expression for the spherical Bessel functions (Equation B.1) into the 


left hand side of the integral above, we immediately obtain: 


1 sin ((v — 15) 








ae t) j,(at) dt = ; B.24 
ff islet) jlat) at = oe (B.24) 
When v — yp, the right hand side has the factor of the form : 
sin ((v — p)= 
im (( #4) BS (B.25) 
CS (Peay ies dee U2 
So we have the following important results : 
ere ; as 
i julat) j(at) dt = Da(Qn +1) (v = p) (B.26) 
= 0 (vy — yt) is even 


x 0 (vy — pt) is odd. 


Using the recurrence relations in the previous section, we can deduce some important 


integrals. Clearly: 

















[FP in(der) in ker) ar = Teen (B.27) 
[inal der) jn-a(kr) dr = Tea (B.28) 
J ings) jnaa(br) dr = TSE (B.29) 
Since 
So Ae easalle)), (B.30) 


kr Qn +1 


therefore, 





[pine tdn) = [> Sa inal) + Jess inna). 
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(B.31) 


Only the j,-1(kr)jn-i(kr) product term will evaluate to a non-zero integral. The 


other term, jn—1(kr)jn4i(kr) will yeild a vanishing integral because n — 1 — (n+1) is 


an even number. Thus: 





















































Oca alae) 1 
nek eS B.32 
I ee OP) > On = Oe) ee) 
Similarly, 
Oe or Ga RT) a 
mater) = B. 
i dr ee InP) = 2e on + Gn +3) ee) 
Using the recurrence relation 
we (nial) = SA (nn  L)incal hr) = ninaa( hr) (B31) 
kr dr S27) = On pT sie eae 
we obtain the important integrals: 
co ay a(n +1) 
n-i(kr) - —— (rjn(kr)) = B. 
i GPP AT) ae an PO) See = yon eT) (Be) 
and, 
Cote ss ae rn 
y dr jn4i(kr) - kee (rjn(kr)) = — (B.36) 


2k(2n + 1)(2n 4+ 3) 
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Finally, using the recurrence relation 


Ta inl Br) = Sag (tinal) — (0 + Deal), (B37) 


we deduce the following integrals: 























[P trincalbr) FS ilk) = ST (B.38) 
and 
co i 8 a(n +1) 
[Oar jnaslkr) Tiny io") = aay Toe a3) (B.39) 











LAs. 


Appendix C 


Associated Legendre Functions 


C.1 Legendre and Associated Legendre Functions 


The @-part of the separated scalar Helmholtz equation has been shown to reduce to 


the following differential equation, known as Legendre’s Equation: 


dP 2x dP 1 m? 
= i.e = 1 
dx? [oetdp hoe (nn ) al? . ) 


The general solution to Legendre’s equation is given as: 








P(cos@) = Cin P”(cos8) + Dinn Qi" (cos) (C.2) 


Here n € {0,1,2,...} and m € {0,+1,£2,...,n}. The functions Q™(x) have loga- 
rithmic singularity at « = +1 i.e., on the z-axis. We will not be concerned with these 
functions since the points on the z-axis, at x = +1, are going to be included in the 
problems. 

There is a certain degree of variation in the definition for the Associated Legendre 
Functions. It is mostly in the definition of the phase that different authors have 
chosen to follow different conventions. In this work, we have throughout followed the 
convention of Magnus and Oberhettinger [22]. This is the same phase convention 
followed by Jackson [22] and the reference handbook by Gradshteyn and Ryzhik. We 


choose to follow the following definition for the Associated Legendre Functions: 








d™ 
Pr(x) = (-1)™(1-27)"P 2 _P,(2), (C.3) 
apt 
(<1) 4 2ymen em 
= — m/2 2 1)" 
>(nly = dxgrtm (2 ) 
where x = cos@ and P,(az) are the Legendre polynomials, which are solutions to 


Legendre’s equation with m=0. 


The Associated Legendre Functions follow these basic properties: 


PoMe) = (DER (C-) 
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Bae) Carte) (C.5) 


Pee). a, 0 (C.6) 
and the orthogonality relation: 


2 (n+ m)! 
In + team eae 





[a2 PM(2) PEC) 


C.2 Recursion Relations of Associated Legendre Functions 


We enumerate a few recursion relations that will be used for solving some of the 


integrals in the next section. Here z = cos# where @ € [0, 7] : 


Pale) — Paa(a) = (n+ 1)v1- a? PPo'(a) (C.8) 
P™ (x4) -—2P™(2) = (n—m41)V1—2°9P™ (2) (C.9) 
oP ge) Sa) =) neton yl = eee) (C.10) 
(1-2) 2 Pre) = (nt 1ePM(2)—(2—mt1)PR(e) (Cl) 
= —nzP?(r)+(n+m)P™,(2) (C12) 

= -V1—22P™"!(¢)— meP™(z) (C.13) 

Cl — 24 pm(z) = (n—m+l1)(n+m)v1l ee eae Cos (C.14) 


dx 
+ mazP"(x) 


(Qn4+leP™(x) = (n—m+1)P%4(2) + (n+ m)P2{(2) (C.15) 


m 
yen 


(Qe) (C.16) 
Equation C.8 can be rewritten by letting n — n — 1: 


Prez) 
VI = 2? 


= —(2n —1) PRG!) + Tale (0.17) 





Similarly, we can write: 


























We can also write: 


=e) Ree 


—(2n —1) P™5\(2) + 


—(2n — 1) Pt y(a) + 


—(2n — 3) Prvg'(a) + 


ae 


E 


rae 
1 — x? 





Pee 


1— x 


SS 





: 


Prog (@ 


1— x 


SS 





; 


PE le 


1-2 


SS 





; 


Pee) 


—(2n +1) P™(«) + 42 





V1 — 2? 


—(2n — 3) Pt a(2) + 3 = 


—(2n — 1) 


—(2n — 3) Prvg*(a) + 


—(2n +1) Pele) He MEL 


Vl—a2P™ (2) = 


V1 — 2? Pia) 


Vilaet Ps) 


Vig? Peg) = 


(2n + 1) 


(2n — 1) 





(2n + 1) 





(2n + 1) 


P™?(a) + 


n-1 





as) 
I 
a 
~ 
& 
iemeaaie 


n-2 
1 — 2? 





i 3 


Bea @ 


1 — 2? 


SS 





! 


Peas 


iE 


1 — 2? 


[Prii(e) — Pha (2)] 


[Pea a(2 


SS 


=P 


SS 
oo 





[Pra (@) — Paa'(@)] 





[Prov(2) — Pay’ (2)] 


£9 


(C.18) 


(C.19) 


(C.20) 


(C.21) 


(C.22) 


(C.23) 


(C.24) 


(C.25) 


(C.26) 


(C.27) 


(C.28) 


(C.29) 


(C.30) 


(C.31) 
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VIF PRCA) = plPt(2) - PAP Ce) (C.32) 

VIRB PB) = Gl) — Pea) (C.33) 

VIF PEMA) = Gy lPtle) — Pala] (0.34) 
Using Equation C.15 we can immediately write the following equations as well: 

2 P™(2) oray (n= m+ 1)PM (2) + (n+ m)PM,(2)] (C35) 
ePEMa) = ey [(m—m + PMN e) + (n+ m— PEGA] (C36) 
PPM) = Ey [(m— myPmi a) + nem + PMe)] (C37) 
2 PM+(2) ory (rm — m + 1PM) + (n+ m + 2)P™1(c)] (C.38) 
2PM) = Gay le—m— DEM) +t mere] — (C39) 
oP) = Gay le—m+ DEM) + (n+ m—2PEF@] (C40) 
rPE@) = Gey le—m + 9PEae) + (nt mPEMa)] (CAL) 

PP) = Geo lm— mPa) +e tm—YEMa] (C2) 
PPC) = ry lin m+ 2PM) + ne m4 EMA] (CAB) 





C.3 Useful Integrals of Associated Legendre Functions 


Now we shall evaluate a number of integrals that use the recursion relations listed 


above and the orthogonality relations of the Associated Legendre functions. These 


integrals are useful for evaluating the scalar products involving L, M and N functions. 


121 


Let us consider integrals of the form: 
T 1 
i dd sind P™*"(cos@) P™(cos) = i dx V1 — @ P™ (7) P™ (a). 
0 -1 
Thus: 
1 1 1 
[eevIqPPROPM) = GIy fe = = ™ P™(2) 


eee 
= “enon to8 
2 


2 (n pee 
(2n — 1)(2n + 1) (n—m)! 











ie dzV1— PPE (er le) = On rE rs 1) “ z = ( m= 0) (C.44) 











[de VIF PEMe)PM2) = ory [ae [P(e - Phala)] PM) 


= orca [de Pa) PM2) 








ie dx JT — 2? P™51(«)P™(a) = oe in 55 . as = (m2>0) (C.45) 











[eve POPE) = Gey fee Pet PRB) — Pate] Pee) 


l 1 
= Gran [Pete Caw cha G) 





, ah bmg (Qn + WED =1) 5 ewe 8) 
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1 | 1 
is de VI — @P™l(2)P™(a) = maple [Pms(e) — Pmi(e)] PR (2) 


1 : 77 in 





[VIM FRB OPE) = ae +3) oe = ( m 2 0) (0.47) 














Integrals of the type: 


= pmtl() P™ (2). 


é 1 
| dO. cos) P™#*"(cosb) P™ ,(cos#) = I, dx Jrwai jen 


We have 


Pe 5 [Ce — myPrte) + (nm — A) Pale): 


Gn =) 
Now since 

PrN) =. on —aypm-te)4 PRO). 

V1 — x? V1 — 2? 
we can continue the recursion on the second term on the right hand side. For every 
recursion the second term on the right hand side will have the ‘n’ index reduced by 2 
where as the ‘m’ index remains the same. Thus eventually we will have the m index 
exceeding the n index and the recursive series will stop at that point. All the terms of 
the right hand side generated by carrying out the recursion will contain factors of the 
type P™ , only. Thus in the product of terms generated by expanding the integrands, 
we have only one term that will be non-zero due to the orthogonality relations of the 


Associated Legendre Functions. 


‘i v m7 m7 
i da Fa Pa) 


= oe dx [(n —m)P™(r)+(n+m— Dee) x 


[-(2n — 3) P™ (2) — Qn—7)P™,—-- | 
(2n —3)(n+m—1 


= fae Peal Pale) 
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(—2)(n +m —1)(n +m — 2)! 





Pri(@)Pea'(e) = 


Les a" tr (Qn—-1) (n—m-—2)! 


(m > 0) 


(C.48) 








ie dx ——— Via Prail (oe) Pea (a) 


= a dx | (n-—m+ Lp ha) +(n+m-— 2)P™3"(«)| x 
n—- =i 


(2 


|—(2n — 3) P52) — (2n—7)PMI = -- 





_(2n—-3 \(n +m — 2) 
— despre a) Peo 
aS ye n—2 n- (x) 


| 








(—2) (n+m-—2)! 

















i - qa 0) Pra (#) = (2n —1)(n-—m-1)! Oe 
[POP 
m ed dx [(n —m)P™(a) +(n+m—1)P™,(a)] x 
[-(2n + 1)P™(«) — (2n-3)P™, ----] 
= a Ore lin —m)(2n + nf (roa ei (ea geal | 
+(n +m — 1)(2n — 3) ‘a dar (a) P(x) 
= (—2) | (n+m)! er 
(2n—1) |(n-—m—1)! (n—m-—2)! 
[de prorat) = (-2) PEPE (m 30) (C50) 











ie 7 ers meat) Pea (2) 
= —_ i a l(n —~m+3)P™5(a) + (n+ m)P™-*(2)| 5 


[-(n — 3) PRG(2) — (Qn -D Pry —-- 


So there is no overlap, and the integral should vanish: 
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[iw Fae) P'(e) = 0 (m2 0) (C.51) 
, d | P™(7) pm 
[ ae n (t)Pagy (2) 
= [ dx P™ (x) [-(2n +1)PMx) — (2n-3)PM, — 
eee nf P™(«) P™(2) 
. 1 a EIR (n+m)! 
ie dx Jira (ey (2) = ra ali (m > 0) (C.52) 
[ew = aOreade) 
= [ dx P™ (x) [-(2n 1)P™ (2) = On =3)P™s = 
[ dx 1 P™(2)P™ (2) =0 (m>0) (C.53) 
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i a Gee ) 


[de BmGMa) [n— 1) PPG a) — On — 5) PRGE- 


—(2n — 1) me aC zt) P" >" (a 


n-1 n-1 











: 1 ing eae oe (n +m — 2)! 
[te POPE) = (DEE (ma0) (Ct 

















[ia + P™(2)P™'(2)=0 (m>0) (C.55) 








Pr (a) Peat (@) 


mo 
rat e 
Qu 
& 
— 
= & 
2 





- L (dt Pa) Pegit(a) — f de VT— a? Po (2) Pas (2) 


(n+m)! i" 2 (n +m +4 2)! 
(n—m)!  (2n+1)(2n+3) (n-—m)!)! 














[ia x P™(a ) ptt) = gi +m)! n+m+2jintm+l)_, (C.56) 


( 
ae ee | Pep eer 9) 





1 2 
fi PP (2) PU (2) 
-1 














126 





: 1 m m+1 2pm m+1 
= - dx Jraat (x)P™ (a) — I, dev l= 2? PO (2) Pe (2) 
a 2 (n+m)! 
(2n + 1)(2n —1)(n—m-— 2)! 











(n+m)! 









































: x? m m+1( >.) — (—2) 
is ie ag ee) Gag (Qn —1)(n—m—2)! a) 
[& F=SePr@Ores') 
1 l ” bed 1 a = 
= ee ds Pa) PEG (2) - ie dx V1 — @P™(2) P™>"(2) 
3 (n +m — 2)! 2 (n+m)! 
(n—m)! (2n + 1)(2n —1)(n—m)! 
1 x? a lg = (n+m— 2)! [(n+m)(n+m-—1) 
de de Fa Pr (2) PIG (2) = 2 | ae 
a de Fs Pa) Pra (2) 
1 1 a és i 1 — bd 
- [a Fao) Pa (@) - [ie Ji — a P™(«) P™'(a) 
——. 2 (n+m)! 
a (2n + 1)(2n + 3) (rn — m)! 
1 x? oe at _ (—2) (n+m)! 
ie 7 loan" (2) Pra (2) = (2n + 1)(2n +3) (n —m)! op?) 











Now we have seen that 


SP (2) = — [nx PM (2) —(n+ m) Pm (x)] 
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and therefore 
T . d 44 
i dé sin @ cos 6 pin (cos 6) | PAY (cos 6) 


1 1 

= nf eas Pee predate j= (n+m) [de Farle oe ) 

+m)! stm \(n+m +1) 
(2n + 1)(2n + 3) 





(n+m-—1)! 
(n—m-—1)! 





(n 
n2 
(n—m)! 


—1] = G+ m2) 





es 


"Spat Ocoee 
| sin cos 6 


ar 


P" (cos 0) PE cos.0) (C.60) 


= 2! 


(n—m)! 





—-m 


n+m)! a +m+2)(n+m-+1) 
(2n + 1)(2n + 3) 











i dé sin @ cos 6 wat (igen 0) P™*"(cos 0) 





# fe FeO RREIN) ~ (nt mn) fdr ea Po) PEN 
(—2) (n+ m)! (—2)(n +m —1)(n+m-—2)! 


= enact een eran) 











ie d@ sin 0 cos @ (Pecos) P™+(cos 8) (C.61) 
0 


a 2(n + 1) (n+m)! 
~ (2n — 1)(2n +1)(n—m — 2)! 














™ . d Vita m 
| dé sin 8 cos @ Ad (cos a) | ag (ce 6) 


= nf du T= P™(2)P™5'(x) —(n +m) ie da Gaara a) Pra (@) 


. el a +m)! 
(2n + 1)(2n + 3) (rn — m)! 
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i d@ sin 8 cos @ (spPeos 0) bie (cos 8) (C.62) 
0 


a (—2)n (n+ m)! 
(2n + 3)(2n + 1) (n — m)! 














. os d TL m7 — 
| d@ sin @ cos 9 (a0 (cos 0) P™>1(cos 8) 


Pe 


n-1 


7 PO(e)PRGM a) — (n+ m) [de as PO (a) PEG) 


n-1 


na qe? 
Tr ———————————— 
is a ee 











_ (n+m— 2)! |[(n+m)(n+m-—1) (—2) (n+m-—2)! 
eer oe | (Qn—1)Qn+1) 1 APO ay (ean 
i dé sin @ cos 6 (5 Psn(co 0) P™>"(cos 8) (C.63) 








_ A(n+m) (n +m — 2)! [in—mtn (ea a) 


(2n-—1) (n-—m!)! (2n + 1) ~ (n+m) 




















ip dx aa nian colar Ges afl dz P(x) |-(2n +1) P™ (2) — (Qn-3)P™, --- | 
1 wh ate _ (n+m)! 
[de GaP PRE a) = (2) PO (C.64) 











1 1 - 
qe dx Vina" (2) Pay (2) 


= ie dePia (2) [-(2n —1)P™51(e) — (Qn — 5) Pmz — -. ] 





1 l 7 Sth he. 
ie dx Vina (2) Pra (2) =0 (C.65) 

















129 


, 1 m m+1 
[ear a ee) 


. [ da P2(x) |-(2n — 3) P™,(#) — Qn —7)Pm,—--| 








a dx —2— P™(2) P™(2) = 0 (C.66) 


Aa af lieegee ® 








[ee SPrloPEV@ 


=f da Ps" ( [-(2n —1)P™7\(2) — (2n —5)P™ - + 


n-1 


= -(2n-1) is _ de PRG) P2G"(2) 





[errors = (C67) 


=i Vl—«2? ” 











Making the substitution n’ = n + 1 we obtain 
1 
a4 de Pa) Pa) 


We have already evaluated such ate aes C.48). Thus we obtain the answer 





toa (2 )Pot (2). 





by simple substitution n — n —1 in the previously found result: 





(—2) (n+m)! 
(2n+1)(n-—m-—1)! (C.68) 





n(a)Py (a) = 


1 x 
d P 
I, A 4 


Similarly, 














[eS Prorrne) 


eed Lake sl i dr 





(2) Pry (2). 


Van 
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Using the result already obtained (EquationC.49) we obtain: 























[ta eee ita) PE Ma) = Ga (0.69) 
iy dé sin @ (ate 0) P™*"(cos 6) 
= =nf dx Vine Pg )Pr(s )-(n+m) [. dz Gaal x) P™*(z), 


The second integral with the substitution n — n+1 is immediately transformed to an 
integral of type given by Equation C.52. The first integral is of the type just obtained 


in Equation C.68. Thus we obtain on substitution: 








[ dé sin @ (3 P™ (cos a) P™* (cos 6) = a : i x ae (C.70) 











A dé sin (ate 0) P™~"(cos 6) 


i if 7 
=nf ds P(e) Pra) - (rtm) f de -—s Pa x) P™"(2). 

Now the second integral with the substitution n — n+ 1 is immediately transformed 

to an integral of type given by Equation C.55 and therefore vanishes. The first integral 


is of the type obtained in Equation C.69. Thus we obtain on substitution: 








(—2)n (n+m-1)! 
(2n+1) (n-—m!)! 





[ dO sin 8 (5 P™(cos 0) P™-1(cos 0) = (C.71) 








: deere) iia 


7, dz P™ (x) |(n-—-m+1)PR (2) + (n+ m)P™ (x)] 


(n nt m) 
= mL dx Peat eee) 
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9 (n+m)! 
3 d7or (ae aa) = (2n —1)(2n +1)(n-—m-—1)! (C.72) 
[. doe (er ila) 
a de PM (2) [(m —m + 1)PMy(2) + (n+ MPM (2) 
= fer oree) 
(n+m-+1)! (C.73) 





7 2 
ie dx xP” DP (es )= (2n + 3)(2n + 1) 


(n—m)! 





a dO sin? 6 (+ P™(cos 0) P™.,(cos 6) 


anf dea Pr (2). P" (2) = (rtm) da Py" (a) Paha (2). 


The second integral vanishes immediately due to the orthogonality of the P(z) 


functions. The first integral is obtained from Equation C.73. Thus on simplification 


we obtain: 





2n 


(n+m-+1)! 








a dé sin? 6 (4 P™ (cos 0) Py (cose) = bee Ones) Ga a 


(C.74) 





I do sin? 9 (5 P™ (cos 0) P™ ,(cos 4) 


anf dra Pra ) P™, (2) — (atm) [ dx P™ (2) P™ (2) 


2 (n+m)! 


2 (n+m-—1)! 








“SNORE Gay She 


(2n—1)(n-—m-1)! 
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‘ a dn oa —2)(n4+1 n+m)! 
[ dé sin’ @ (a0 (cos a) P™ (cos @) = a 7 wos 2 i a ~ a (C.75) 














C.4 Additional Integrals 


The (Associated) Legendre’s differential equation can be written in the following form: 


2 


d Boao m py: 
zla-2) Z| + [untae PO, (C.76) 








where P” are the usual Associated Legendre functions. We can write 





























d ae (2 p™ dp” 
— (1 — 7? ees Soe no_9 n 
a 2 aa are de? 
and 
é dP d ae dpm dp 
ora i 2 be iss — eae {= 2 n —— 2 n n 
So 
dP | _ d dee dP™ dpm 
Pr — 1 2 7 = —_— | 2 {pur me — (Po 2 n ni : ; 
2 | a) “| App | le erreurs (C.77) 


Multiplying both sides of Equation C.76 by P” and P” respectively, we obtain: 


2 




















d dPp™ m 
Pe = 2 le = PEP = 
we [a2 SE] + fain 41) - | pres =o 
d dP” m? 
m = 2 n i i ‘le = mM m 
pee [a 2) 2] + [wor +1) - | ree =o 
Or, 
d dPp™ dP™ dP™ m? 
“Ja —22)pr@2| ~a -_2)42 1)- —" | pm pm — 
dz ( w)Pr =| ( 2 dx dz + [ge + ) a ade 


(C.78) 

dP? dP™ m? 
= _ pt ni n / / =: mpm _ 
| (1 «\—, ce + [rent 42) | eee 0 


(C.79) 


UL, 
dz 











& |a-2ye 
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Subtracting Equation C.79 from Equation C.78 we obtain 








d dPp™ dP” 
[a2 {pe S2 — pe SEN 4 n(n 1) nl 4) Pre = 0 


Integrating both sides with respect to x, between -1 and 1, we obtain 
Oo ae 
1 —2?)¢ pu —* — pm — 

= 2) | pe SE — pe | 


The first term vanishes at « = +1 since the factor (1 — c") vanishes. The second term 


1 





1 
+ [n(n+1)—n'(n' + Wf iio eke peal 
-1 





-1 


must therefore be zero. So, ifn # n’ then we arrive at our usual orthogonal relation: 





[de PR(a)PEa)=0 (nd rn) (C.80) 











Adding Equation C.79 and Equation C.78 we obtain 





a2) pe SE + pe SEN 4 tnt) tl + Les 





dx dz dx 
dP™ dP™ 2m? 
( are dz tilt 


Integrating both sides with respect to x, between -1 and 1, the first term on the left 
hand side will vanish again because of the (1 — 2) factor. Therefore we immediately 


arrive at the following relation: 


Lt n(W’ $1) p dpm dP? 
mene Daf aw prem = | ae {== aban ae preg}. 
-1 -1 








2 io dx dz a ge ee 


Changing the variables to @, so that 





d Lod 
— p™ == — P™ (cos 6 d 
dei *) moe oo ae 
d d d d 
(1 - w*) Pe (2) 5 Pe (2) = aoin (08 8) agin (cos 6), 


we arrive at the important integral 


wv 2 
| dé sin @ (sree 0) P(cos 6) + 5 Pi (cos 6) Pr (cos 0) 
n(n + 1) + n'( 

2 





PED fae PR(e)PB(2) 
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Using the orthogonality of the P’, this immediately implies: 











n ; d d m? 2n(n+1)(n+m)! 
0 6 {| —Pp™— Pp™ m pm | — er, C.81 
j ace (5 fr do" sin” gi" fr 2n+1 (n-—m!)! ( ) 
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Appendix D 


Evaluating Inner Product Integrals 


D.1 Introduction 


In this appendix we show the details of carrying out various integrations, required to 
obtain the inner products in the determination of the expansion coefficients of a plane 
wave. We will employ various results summarized in some of the other appendices, in 
particular the appendices on spherical Bessel functions and the Associated Legendre 
functions will be used frequently. We have: 
a SH (2s + 1)j,(kr) 1% (ove) P!(cos a)e"8 P!(cos Ae"? 
(s+ 1)! 


s=0 1=0 








(s— I)! l 18 pl -il 
+ P'(cos ae"? P!(cos 6\e7"!? ; D.1 
yo SFPileos ae (cos) (D.1 
where a, / specify the direction cosines of k, and 9, ¢ give the direction cosines of r. 
Also |k| = & and |r| = r in the above expression. The unit vectors of the cartesian 


coordinate system assume the following form in spherical polar coordinates: 


e, = sin@cos¢e, + cos4cos deg — sin dey (D.2) 
e, = sin@sin ge, + cos@sin deg + cos Gey (D3) 
e. = cosée, —sinfeg. (D.4) 


The e** term involves a double infinite summation in the indices s and J. As we 
shall soon find out, the various orthogonal properties of the functions involved in 
these integrations of the inner products will eventually allow only a few surviving 
terms. This will allow us to carry out these integrations analytically, giving us closed 


form expressions for the final results. 
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D.2 Evaluation of (E,|Mj nn) 


These inner products are necessary to evaluate the a%,, coefficients. We recall the 


Th 


definition for the inner product: 
fore) T 20 ‘ ; 
(E.|Mynn) =i ar [ dosind | dé e, -M*, ei =), e, -M*e**,  (D.5) 
0 0 0 rod 


where the r@¢ below the last integral sign is a short-hand notation indicating that the 
integration is to be performed with respect to all three variables r, 9 and ¢. Similarly 
if we had only r@ below the integral sign, it would indicate that the ¢-integration has 


already been carried out. Now, 


Pe 0) _. dp™ 6). 
M* = —imj,d kr) LED) --intg, 2 jn kr) Pees) nd. (D.6) 
The M*,,, does not have any e, component, therefore we obtain 


(E.|Minn) = : (- sin 6 eg) . egg 


et eg) ete 
10 sin 6 


im | kts (kr) P™ (cos 6)e*. (D.7) 
rod 


This clearly vanishes when m = 0 => az, = 0. The e*** factor has ¢-dependence of 


the form e+", Therefore, the ¢-integrations involved in this case are: 
20 ; : 
i doee™? = QW bm and, (D.8) 
0 
20 : ‘ 
| dé ge ee = Oe a ene (D.9) 
0 


Thus among all the terms in the /-summations, only / = +m terms survive as a result 
of the orthogonality properties of the complex exponential functions. Since / > 0 
in the summations for e"**, when m > 0, only the terms corresponding | = m will 
contribute, that is those terms containing the factor e~"”’. On the other hand, when 
m <0, only the e’”8 terms can be non-zero. The 6-integrations involved are of the 


form : 
i dé sin 6 P!(cos 6) P™ (cos @) (b= cera | (D.10) 
0 


Until we know the constraints imposed on the values of s, we still have another 
potentially infinite summation to deal with. Let us consider the r-integration for the 


moment: 


o dr j(kr)jn(kr) =0 (s—nis even, sn) (11) 
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Let us consider the parity of the @-integral. Since P™(—x) = (—1)™*" P(x), we 
immediately see that the parity of the integrand in the theta integration is: 


1 
/ dx P\(z)P™(x) has parity ~ (—1)st"*™, 
ai 


Since / = +m the 6-integral will be non-zero only for s +n equal to an even number, 
ie. +n = 2k, where k is an integer. Equivalently stated, the 9-integral may not 
vanish only when s and n differ by an even number, including zero. On the other 
hand, the r-integral vanishes when s and n differ by an even number, except zero. 
Thus the only possible term that survive from considering the parity of both the @ 
and the r integrations is the n = s term. Thus, from the double infinite summation 
we have only one surviving term that can contribute to the integral for the inner 
product. Thus with /= m and s = n, the @-integral is evaluated in the following way: 

2 (n+m)!)! 
2n+1(n—m)! 





Tw 1 
| d@ sin 6 P™ (cos 6) P'™(cos 6) = i dx P™(2)P™(«) = (D.12) 
0 -1 


and the r-integral becomes: 
i dr ju(kr)jn(kr) = 
0 


Thus, for m > 0, we have: 


T 


2k(2n + 1) 
(E.|Minn) = im | ekt 7, (kr) P™ (cos 6)e"? 
rod 
foe) ™ 20 On 
2 im | ar | dO sind f dé So i#(2s + 1)j,(kr) x 
0 0 0 a) 


5 — /])! ; ; 
2 ° i m P!(cos ae“? P!(cos O)el? 
=o" 








+e zs 5 P, Pi(cosae" Poste} i hr) PE (co bye? 
(s 


= Dri (2s +1 oe dr j,(kr) in(kr)) x 


= ! . “ 
(=m PD (cos ae" [” a8 sin 8P2 (cos 8) P2™(cos 0) 
(s +m)! 0 


= 2mm 2"(2n 41) (Since only s = n allowed) 


= x 
2k(2n + 1) 


Sm) P" (cos ae? 


(n+m)! ” 


2 (n+m)! 
(2n +1)(n-—m)! 
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Thus the expression for the inner product becomes: 





mr? 2m 


ee eon) 


P"(cos aye? (D.13) 











D.3 Evaluation of (E,|N nn) 


The (E,|Nimn) integrals are used to evaluate the b%,,, coefficients. We have: 


Th 


fore) wT 20 : : 
(E.|Nmn) = | dr | dO sin 0 | dé e, -N*,,e** = | e,-N*e**, (D.14) 
0 0 0 rod 

















Now, 
In( ke 1 atk Ke 
Nr, = n(nt+ {V2 *) pe (cos bee, + i Arsak r) dr {cos ee, 
1 d(rjn(kr) P™(cos@) _,; é 
“ a oe byt 
Ee ude sind “? (Pa) 
Therefore, 
(hk . 
e,-N*,, = n(nt 1yrat t) cos 6 P™ (cos 0) e*”"? 
r 
1 d(rjn(kr) . ,dP™(cos@) _,; 
— g— wae. Da 
kr drtCSOC eae) 
As in the previous case, the ¢-integrals will be non-vanishing for 1! = +m only. 


However, the r-integrals are now of the form: 
Led(rjaCh 
fe dr j5( kr) a ye dr j,(kr) — mate ( dae), 


dr 
As derived in the appendix on the spherical Bessel functions, such integrals are non- 





vanishing when s =n +1 or when s — n is an even number. Let us now consider the 


parity of the 6-integrals: 
i dO sin 8 cos 6 P™(cos@) P™(cos@) ~ (-1)it™trtmts J (_y)strtt 
0 


[ dé sin @ sin 8 (Peoos 0) P™(cos 6) ~ (aqimretitets tas (aye, 
0 
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where the sin@ does not change sign (since sin@ = sin(a — @)) and the parity of 
the derivative of Associated Legendre functions will be one greater than the corre- 
sponding Associated Legendre function (see eqn C.11). Thus in either case above, 
the 9-integrals will vanish when s + n is an even number (or equivalently, s — n is 
an even number). The 6-integrals may be non-vanishing when s = n+1. Thus, 
by considering the mutual parity of the r and @ integrals we conclude that only the 
s=n+1 and /= +m terms can contribute to the expression for (E,|Nmn). 


For m > 0 we have the s = n+ 1 terms: 


272 +4 (2n + 3)e LL eeray greron x 
jn(kr) 


(m4 +1) Le dr jn4i(kr) ie ) ([ dé sin cos@ P(cos 8) P™ (cos 0)) 
0° 1 d(rjn(k ™ d 
= (/ dr jnaa(hn) ge SERED | . (/ dé sin? 6 qin (08 6) P™ ,(cos a) . 


The 6-integrals are readily obtained from Equations C.73 and C.74. The r-integrals 
are obtained from Equations B.33 and B.36. Thus we obtain the s = n+ 1 term: 





(n+1—m!)! 
(n+1+m!)! 


(14 + Rone y(n i 5) (= a in 3) i) 


- (- 2k(2n =r + 5) , (a + Hon + 3) al 


ee a ae 
nt? 2n(n—m-+1) enim pm 
k (2n + 1)(2n 4+ 3) 


Qn i"t (In + 3)e"? P™ 


bi (COS a) 











= grt 





(cos @). (D.17) 


Similarly, the s = n — 1 terms are given by: 


(n—1—m!)! 


(n-—1+m)! ~ 


Ck dé sin @ cos @ P’ (cos @) P™ ,(cos 0)) 
0 


Qn int (2n — Lje*™? pe 


(cos a) 


(m4 +1) _ ar ju-a( bre 


= ts ir jnaalhr)g CRED ) : ([ dé sin? 6 P™(cos 8) PE (cos a) 


Now the 6@-integrals are given by Equations C.72 and C.75. The r-integrals are ob- 





tained from Equations B.32 and B.35. Therefore we obtain for the s = n —1 term 
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: =) =) 
Qn i”! (2n — Le"? P™ , (cos ae x 


(4 + 1) en = Te za 5) (os; = in +1) a eee 7) 
7 (- oie, + 5) ) (-a STTEST (n ea 7) 


= oot 2(n as 1)(n a m) e7imB pm 


k (2n —1)(2n +1) ee 


The sum of the s = n+ 1 and s =n-—1 terms give us the desired expression for the 

















(cos a). (D.18) 


inner product (E.|Nimn): 





Qn? etme n(n —m-+1) pm 
k (2n+1) Qn4+3 "tH 
_ (n+1)(n+m) 
2n — 1 


(E.|Ninn) = bee 


(cos a) 





(D.19) 





P™ (cos a) 











D.4 Evaluation of (E,|L,,) 


The (E,|L,,,) integrals are used to evaluate the c?,, coefficients. These can be eval- 


uated in a way very similar to the process of obtaining the (E.|Nimn) integral. Now, 
dj,( kr) In(kr) dP™ (cos yy. 
d(kr) kr dé 


-—imd 








Lyn cos @ P™ (cos @)e7'"? — sin 0 
(D.20) 


Again, by similar considerations on the parity of the ¢, 6 and r integrals, it is easy to 
see that only the /= +m and s =n+1 terms can be non-zero. For m > 0 we have 


the s=n-+1 terms: 


1, -—im m (n aT We m)! 
kQrt +t (2n + 3)e ech ars grees x 


([° dr ial gen) (| dO sin 8 cos 0 P™(cos 6) P™ ,(cos 0)) 
-({° dt jnga(kr)? mir) ): (/ dO sin * P= (cos 0) Ps (00%8)) 
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The 6-integrals are readily obtained from Equations C.73 and C.74. The r-integrals 
are obtained from Equations B.33 and B.39. Thus we obtain the s =n+1 term: 


n -im m (n ar iS m)! 
kQre A (2n + 3)e eC) pare x 


(Grease ym (cesvceemce as 


(sca =f, nes + 5) , (or tg ne +3) we nie ) 


(n M+ 1) e7imB pm 
(2n + 1)(2n + 3) Eas 











= Ir? grt 





(cos a). (D.21) 
Similarly, the s = n — 1 terms are given by: 


| Suite 
kQr gee (2n _ jen" P™ | (cos pee a 
(n-—1+m!)! 


( i ae inal gpl) . ( | " dO sin cos @ P™ (cos 6) P™ (cos 6)) 
-(f° dr jna(kr)2 ud ): ue Aisin 9 SPM (cos a)Pr (cos). 


The @-integrals are given by Equations C.72 and C.75. The r-integrals are obtained 
from Equations B.32 and B.38. Thus we obtain the s = n—1 term: 





2k(2n —1)(Qn+1)} \(Qn —1)\(2n +1) (n—m-—1)! 
(serjerrn) (a as 


+1) 
» (sea = ne s 5) (- (2n “ oe 1) = dl 


(n —m) e-imp pm 
(Qn —1)(2Qn +1) “Fr 


k2Qni™1 (2n — le? Pp (cosa) 














= Qn 7 m—1 


(cos a). (D.22) 


The sum of the s = n+ 1 and s = n-—1 terms give us the desired expression for the 


inner product (E.|Linn): 
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-—imp 
(Es\Le,) ea" 2a" a) P™ (cos a) 
(2n+1) | 2n-1 


(n—m-+1) 
2n +3 





(D.23) 


Pr (cos a) 











D.5 Evaluating (E,|Nmn) 


We write e, in the following form: 
e, = ; isin 0 (ef# + e7'#) e, + cos 9 (ef# + en) ep +2 (ci# = eWi#) e,| . (D.24) 


Thus we have: 








7 1 jn(kr) . m id -id\ .-imd 
Nise = 5 [eee = sin 6 P™ (cos 6) (c +e Je 
ll d(rjn(kr) dP (cos 9) id -id\ ,-imd 
gee a if Je 
1 d(rjn(kr) P™(cos 8) 7 ig -ia\ .-imd 
nr a a a wm : D.2 
ep dr sin 6 (¢ a © (29) 


The ¢-integrals assume the following form now: 
20 : : 5 : 
ij dé (e'@ + ew) embell6 _ In (Comet? Shige) (D.26) 
0 
20 : ‘ : : 
i dé (c’# + e¥#) ete te 2 Oe (Ofna Oronea) 3 (D.27) 
0 


Thus only 7 = m+1 terms are allowed from the ‘summation in e***. The 6-integrals 


and their respective parities are as follows: 


a dé sin 6 P” (cos 6) sin 6 P™*1(cos 6) ~ a aaa a (—1)"ts#1 
0 


as 7 g 
| dé sin @ cos@ as 0) 3s Calera a (Sian 
0 
(cos 6) m1 m+n+tm+l1+s n+st1 
‘a no) 2 Pm '(cos8) ~ (=1) (—1)nts#1, 


Thus only when n+ s +1 is even can these 9-integrals assume non-zero values. The 


r-integrals are of the form: 


es arpa (kr)? 





a(k 
a fein = a): 
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Again, these r-integrals can be non-zero only when s — n is even or when s =n +1. 
In the former case, the corresponding 9-integrals will be zero. Thus the only allowed 
values in the double infinite summation are terms corresponding to 1 = m+ 1 and 


s=ntl.ia 


Now the s = rn+1 contribution is given by: 


Qn i"*(Qn + 3) rn 451) (a dr ineathry) x 


1— 1 a x 
: = — + i eW*m—V98 p™1 (cos ax) : dé sin 8 P™ (cos 8) sin 6 P45" (cos ) 


(n+1—m-1)! 
(n+1+m+1)}! 


+ ( [te jnthr ge EOD) 











este pri (cos a) ‘ie dé sin @ P™ (cos 8) sin 8 P™4" (cos a) 











(eee eo UP pe (cosa re do sing EN”) Flees) cos 6 Pas (cos @) 
n m— 
—~m—l1)!. 7 
7 : m 7 rye tn at oa a) | dé sin 6 GE NCO) og mA oi a) 
n m : 7 
ae 1 d(rj,(kr 
+m (/ dr inea(hr) ge CED ) x 
(nt+1—m+1)! es cose oy 
(in a ie ee ee ae vom ae (cos a) Le a sing n@ : Pr (cos 6) 


(n+1—m-1)! 


ee oe) 
(n+ltm+l)- 


—i(m+1)6 endl fe m+1 
Pry (cosa) | dé sin @ arr ee Pry (cos ah : 





Using Equations C.45, C.47, C.60, C.62, C.64, C.65 for the 4-integrals and Equations 
B.33, B.36 for the r-integrals, we can evaluate the s=n+1 term: 





ri?t(2n + 3) a +1) (san 7 Tyan : 5) x 








(2n—m+2)! yg nm-1 2 (n+m)! 
(n+m)! x OP (08a) (os + 1)(2n + 3) (n—- =) 
(n—m)! _, eaaeaeat 2 (n +m +4 2)! 
"Ga eae OP ae (co) Gersvees) (n—m)! )} 





_ (- 2k(2Qn aE + 5) 
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joa =m + 2) -iim-)8 pm-t(cos a) (- Teg et | 
(n+m)! (2n + 1)(2n + 3) (n — m)! 
_ ! : 
pm ETS Pa (cosa) x 
ntm ! 





fesat(Segriaeese-»)} 





7 (- 2k(2Qn on + 5) je 


- Dh ag 
{enn roo EM 
n my). 


: aoe eo (™+1)8 P™1 (cos q) (-er) \. 


On simplification, this gives us for the s = n+ 1 term: 


2 
gots nr 


k (2n + 1)(2n + 3) 





[(n —m+2)(n —m + 1)P™51(cos ae ("VP 
— PRE (cos ae mtv? | 

Similarly, the s = n — 1 term is: 

Qn i1(2n — 1); an +1) (s dr inal E) x 


kr 
(n-—1—m+1)! 
(n—1+m-1)! 
—l—m-1)! 
ge m ) 
(n—1l+m+1)! 


ip ts ir jnaalbr)g CRED ) x 








e~(m-1)8 P™=1 (cos a) [ dO sin @ P™(cos 6) sin 8 P™>"(cos 8) 


n-1 





n-1 


eW(m+1)8 pm+1 (cog a) [ dé sin 8 P™ (cos 6) sin 6 P™41(cos a) 


kr dr 





eee le : ™(cos 6 
{te m + UY! -itm—1)8 pm=1( cos a2) [ao sin TCS) cos Pm! (cos 8) 
0 








(n—1+m-1)! 
—-l-m-1)! _, ™ dP™ (cos 6 
ale * = eo — zt i eW(m+1)6 pm+1 (cog a) | dé sin @ ah(cos) cos 6 P™+!(cos a) 
oo ; 1 d(rjn(kr) 
+m (/ dr jn-a(kr 7 x 
(n = dat Ly —i(m-1)8 pm-1 “ : P" (cos 9) m—-1 
te Sarat)! € P™>*(cos a) [ dé sin ae P75" (cos 9) 
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ee) ee ™ (cos 6 
= fu i ) ee prt) (cos a) f d@ sin 8 Ex (eos @) 
(n—1l+m+1)! i 0 sin 8 





Preiveos6)} | 


Using Equations C.44, C.46, C.61, C.63, C.66, C.67 for the 9-integrals and Equations 


B.32, B.35 for the r-integrals, we can evaluate the s =n —1 term: 


ri" (Qn — 1) rn +1) es - ek : 5) x 











(n—m)! e-i(™-1)8 pm-1 (cog 9) (— 2 (n+m)! 
ieoh Frat ) ( (2n —1)(2n + 1) = 
(n—m—2)! _, aaa as 2 (n+m)! 
+ (n+ m)! en Pav (cosa) (opm) 





% (ste + 5) ‘ 


2 = my —i(m— ma 
Nate 2 m9? PG *(cos a) x 


(Matai Cnt Chai Sea) 


(rn —m — 2)! eae aaa 2(n + 1) (n+ m)! 
(n+ m)! en ce (cosa) (ata 














is (sean a3 5) | 
Lice ae a) x (es) 
n—m—2)!} 


Gam ee a) } . 


On simplification, this gives us for the s = n — 1 term: 


Suet (werd) 
"  & Q@n—-DQn4)) 





[(n +m)(rn+m— 1)P'™>1(cos eye en? 
— P™*(cos aetmtn8) 


Finally, adding the s = n — 1 and s = n+ 1 contributions we obtain: 
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gut 1? 
E,.|N ——__ — 
alee = ea aay a 
n(n—m+2)(n—m+1) ates 
Gee Neale 
nr m -—i(m 
“oregon (D.28) 
Nia btm =D peng 
Tr — 
(n = 1 m -—i(m 
~ PED prt cos ayer? 











D.6 Evaluating (E,|M,n,) 


We have 
* = ohne P" (cos 8) id -id\ .-imd 
M*,°@: = 5 —im jn(kr) cos @ era as (c +e ) € 
55 dP™(cos@) ¢ ip -i6\ .-imé 
—t jn(kr) =a (c —e ) € . (D.29) 


As in the previous case, the ¢-integration allows only the 1 = m+1 terms to be 
non-zero. The 6-integrals in this case are of the following forms, indicated along with 
their parities: 

a P™ (cos 6) 


: d@ sin @ sin 


T ™m g 
j d6 sin g LE (608 9) pomsit(egg 9) ~ (—1)mbntibmattts . (_yysts 
0 


cos 6 P™*!(cos6) ~ (—1)ymtntitmeits ~ (—-1)"t* 


So, only the s = n and s +n = even terms are allowed by the @-integration. The 


r-integration in this case is the straightforward integral: 


Ss dr j,(kr) jn(kr), 


which is non-zero only for s = n or s +n = odd. Thus, the only allowed term from 
the s-summation when both r and 6-integrations are considered together is the s = n 


term. So for m > 0 we have (Ez|Minn) equal to: 
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1 
QrittQn+1 )5 sf dr jn (kr) jn(kr)) x 


— 1 
mB mt iy i(m—1)8 p™—1 (cos a) A dé sin 6 Histos@) cos 6 P™~*(cos ) 
(n+m-—1)! sin 
— ee (m+1)B pm™+1 (cos yy dé sin @ Pecos) cos 6 P™* (cos 6) 
n+m 
_ 6 
28 Se), Sim ae ae ‘(cos a) if dé sin @ dF (C08 #) pm—1 (65 6) 
(n+m— 1 dé 
Le aN 
(rte), ponerse Ie peer y (cos a) ve d@ sin @ dF (COS 9) pm (695 6)|. 
(n+m-+1)! dé 


Using Equations C.68, C.69, C.70, C.71, for the 6-integrals and Equation B.27 for 


the r-integral, we obtain : 


(Bi|Minn) = 2"t1(2n +1) (ace) x 








nes (n—m-+1)! en iln-1)8 pm=1(cgg a) _ An +m—1)(n+m-—2)! 
(n+m-—1)! (2n + 1) (n—m)! 

_ mn —m— LY -ilm4)8 pm+1 (cos % _2(n+m)(n+m—1)! 

(n+m +1)! (2n +1) (n—m-—1)! 


n—-m+1)! 


_( 1)! 
~(n+m—1)!- 
a eae 
"ae iD 





en iln—1)8 pm eee Qn een) 


(2n+1) (n-—m!)! 


2m core tone (Bett ot 





On simplification, we finally obtain: 





gntl 1? 


(El Mia: oS Gn 1) & [Pe (cos jenni 
n 


+(n+m)(n—m-+1) Pus €6s aetm—1)8)] 


(D.30) 
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D.7 Evaluating (E,|L,,,) 


We have 


1 
Linn @s = ks (agit) sin 6 P”"(cos 6) (ci@ + e~i#) ee 








jn(kr dP” (cos 8) ib -~id\ ,-ime 
+ ie cos 9 78 (c +e ) € 
In( kr) Pi (cos 6) GN Aas 
—e' ic D.31 
a kr sin 9 (e (D3 ) 


As in the case of (Ez|Nmn), the ¢-integration allows only the 1 = m+ 1 terms to 
be non-zero. The 6-integrals in this case have the same form as their conterparts in 
the evaluation of (E,|Nimn), so that only when n + s + 1 is even can these terms be 


non-zero. The r-integrals are of the form: 


ie dr j5( (kr) 2 


and these integrals are non-zero for s = n+1 or s+n+1 is an odd number. So in 





a fe dr j (kr) ain 


this case too, the combined action of the # and r integrations allow only the = m+1 
and s = n+1 terms to survive. 


Now the s = n+1 contribution is given by: 


on it1(2n + 3)= ala ir inaal br) jul) x 


(n+1—m+1)! 
ee 
(n+1—m-1)! 
(n+1+m+1)! 


(feat), 
(n+l1—m+1)! 

eer ‘ 
(n+1—m-1)! 
(n+1+m+1)! 


+m i“ ar inaalbny”) x 





er lm )e pme1( cds a) im dé sin 0 P™(cos 6) sin 8 P™> "(cos 6) 





errr prt (cos a) [ dé sin @ P™ (cos @) sin 6 P™4" (cos a) 





P™ (cos 6) 





eC sae ce gi (cos a) fr dé ep tn (0088) cos 6 P'™> "(cos 6) 


. ™ dP 6 
enrtetl)e prt (eos) | dé sin 8 ars(cos) cos @ P™41(cos a) 
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jee ge lm-1e La elkare cos a@ if dé sin pe cOse) = pet ‘(cos 6) 





(n+1+m-1)! . ne a 
ii d 
eae Geto -i(m+1)8 p™+1 (cos cy) fe Ae ain 8 (cos) Pritt (cos 8) >|. 
(nt1+m+1) sin @ e 


Using Equations C.45, C.47, C.60, C.62, C.64, C.65 for the 6-integrals and Equations 
B.33, B.39 for the r-integrals, we can evaluate the s =n+1 term: 
+1) 
ka iPth(o ee m(n 
an n+3){( Dk(2n + 1)\(Qn +3)) * 
(rn —m + 2)! Bae eee 2 (n+m)! 
(n+m)! : Pon (cone) (2n + 1)(2n + 3) (rn — m)! 


(n—m)! 2 (n+m +4 2)! 
Geel 8) ( === (n—m)! )} 








-—i(m41 m+1 
€ ( be gies ( 








=: (ae +DQn+ 5) ‘ 
(n —m + 2)! Pe eee 2n (n +m)! 
(n+m)! Prva ( ) ( (2n + 1)(2n + 3) a 
(n—m)! 


a Gt eae 





eer? pe eos a) Xx 


(= +m)! (“* +mt+2)(n+m-+1) _ m)) 


(n—m)! (2n + 1)(2n + 3) 








a (se Ga + 5) * 
{eon ee ~i(m—1)8 pm=1 (cos a) x (0) 


eee (n—m)! seta cont 2(n + m)! 
cer Me acer il 


On simplification, this gives us for the s = n + 1 term: 
2 
n+ aq pmti -i( 
a (Qn + 1)(2n +3) Fpl (COs ae 
—(n—m+2)(n—m+ 1) P5"(cos ee ten —8) ; 


m+1)8 





Similarly, the s = n — 1 term is: 


Qn i"-1(2n — 1)— all f dr jn1(kr) iim isk) x 
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et sh ee 7 
te inissl) eg Tn se pet (eds a) [ dé sin 6 P” (cos 8) sin 6 P™>\(cos 6) 








(n—1l+m-1)! 
(n-—1—m-—1)! 25 +1)6 41 wT . ; ii 
esr Pe (GOs a) | dé sin 6 P”(cos @) sin @ P™*"(cos 9) 
(n—1l+m+1)! aod 0 Ms er 
ii (a ir jnaalbr =) x 
0 kr 
te a= — + 7 en im—1)8 pm-1 (cos a) | d@ sin @ Geet) {cos cos 6 P™7'(cos 8) 


vee) 
(n—1l+m+1)! 


+m (s ar inaalbry) x 


eT et 1) ie me ‘i 2 PE (GOSO) 2% 
te =p)! € P™>*(cos a) | d@ sin 6 Sarr ee P™>51(cos 6) 
oe ee a hee 7 m 
me ma) e~(m+1)8 pmt1 (cos a) | dé sin 6 den (COE) ico 9) 
(n-1+m+1)! 0 sin 6 


m ™ (cos 6 
e~Km+1)8 p™*1 (cos a) i, d@ sin @ aces) 
0 





cos 8 P™*!(cos a) 











Preiveos6) | 


Using Equations C.44, C.46, C.61, C.63, C.66, C.67 for the @-integrals and Equations 


B.32, B.38 for the r-integrals, we can evaluate the s = n —1 term: 











Bene) (seas ah + 5) * 
(eazy or ©) (-g ype mm) 
ee nmi (cosa) (a 7 Gn +1) = a) 





= (saan = ne + 5) ‘ 


(Ram)! sees on 
eon 1)8 P™=1(cos a) x 


(ae (n-m+n ("2 -*=))) 


(2n —1)(n—m!)! 2n+1 n+m 


(n-—m— 2)! -itm-+1)s m1 (cos o 2(n + 1) (n+m)! 
= (n+ m)! Pry ( ) (pt 
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2 [sa =H ri 5) ‘. 
coe ae ene opal cos a) x (- 
5 (n —m— 2)! 

(n+m)! 


2(n+m— =I) 


(n—m)! 


e~lm+1)8 pm+1 (cos a) o} é 


On simplification, this gives us for the s = n — 1 term: 
nr 

(2n — 1)(2n + 1) 
+(n+m)(n+m— 1) P75" (cos ae ten? 


“n+1 _ pmtt —i(m+1)8 


mt (cos ae 





Finally, adding the s =n —1 and s =n+41 contributions we obtain: 

















grtlz2 
(2n — 1)(2n + 1)(2n + 3) 
—(2n —1)(n-—m+2)(n—m-+1) Pi (cos ee ere 
+ (2n —1) pms" (cos are (m1)8 (D.32) 
+ (2n + 3)(n + m)(n +m-—1) P™> ‘(cos ae tm—1)8 
— (2n +3) PP" (eas eje mre 
D.8 Evaluating (E,|Nmn), (Ey|Mmn) and (E,|Lmn) 
We write e, in the following form: 
e, = sin§@ singe, + cos@ sin deg + cos deg (D.33) 


= > isin 0 (ci* — a) e, + cos 9 fee - e~i#) €g +2 (ci@ + e“#) e5| 
When compared with e, which has the form 
e, = ; isin 0 (ei# + ei#) e, + cos 9 (cf# + er) ep +2 (ci# = ei?) e,| , (D.34) 


we can immediately notice that the only difference between the expressions for e, and 
e, is that the former can be obtained from latter by changing the sign of the e~*? 


terms and finally multiplying the entire expression with an additional factor of —2. 
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When we observe the ¢-integrals: 


20 “ : : : 
| dé (ci* + ene) eee? 2 OR (61m—1 £ é1m41) (D.35) 
0 


20 7 ‘ : : 
i. dé (c# + ene ete = Or (Ohegaa ae Open 2n)9 (D.36) 
0 


it is clear that the above transformation will lead to a change in the sign of the m+1 


term, while the m — 1 terms will not change their sign. The only difference between 


the E, and E, terms arise from this source, since the other factors in the integrands of 


(Ey |Ninn); (Ey/Minn) and (E,|L,,,) are the same as the corresponding E, integrals. 


Thus we can obtain the E, integrals from the corresponding E, integrals by changing 


the signs of the ‘m+ 1’ terms and finally adding an extra factor of —:. Thus we obtain 


the following expression: 








(Ey[Nmn} = 7 a\q * 


zn Wr? 


(2n +1) k 
n(n-—m+2)(n-—m41) oe 
7 a ) pm-1(cos axe (m—1)8 
n 2 ms 
+ Bap ay Patt (cos ale (m+1)6 oa 
(2n — 1) 
(n +1) 


+ (on 1) et (008 ae nr? 
nr omemst 








teeing (cos ees ne 











(Ex|Lmn) 





aa? 


Qn —1)Qn+1)Qn+3)~ 
—(2n —1)(n-—m+2)(n—m-+1) Pe (ees ae t™—1)8 


— (2n — 1) P™4*(cos axjetm+1)8 (D.38) 


+ (2n + 3)(n + m)(n +m-—1) P™>1(cos ae tm—08 
+ (2n +3) Po" (eos axe tm+1)6 


n-1 
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2” Wr? 


(Ey|Minn) = Ona) &. I- Brtt(cos ae fmt 


+(n+m)(n—m-+1) P”~'(cos ae tm) 


(D.39) 











D.9 Evaluating (Minn|Minn), (Nimn|Nmn) and (Lmn|LEmn) 


These integrals are easy to evaluate by using the expressions already derived in 
Chapter 3. Since we are interested in a plane wave expansion that is regular at the 
origin, i.e. an expansion in terms of the spherical Bessel functions alone, z,,(kr) repre- 
sents j,(kr). The 6 and ¢ integrations are already performed, and the r-integrations 
can be easily performed using the integrals derived in Appendix B. Since we know 
from the orthogonality relation between the M functions (Equation 3.22) that 

7 A4rn(n+1)(n+m!)! 


20 wv 
d i d0sin9@ Mmm, -M*,,, = Gypsies “Ce 
[ ? 0 yt Ben 2n+1 (n= myn! r) cDa0) 





therefore using Equation B.27 we can integrate both sides of this equation giving us 


(with n =n’ and m =m’): 





2x? n(n +1) (n +m)! 
(Minn| Minn) = k (2n+1)? (n—m)! 





(D.A1) 











Similarly, using the orthogonality relation between the N functions, (Equation 3.24) 


where nown =n’ andm=m/’': 


20 wT 
| dé | d0sin® Non ~Ney1 
0 0 


4rn(n + 1)(n +m)! 
(2n+1)? (n—m)!)! 





[(m + 1)52_1 (hr) + 2324,(kr)| bnnSmm', (D.42) 


and the integrals in Equation B.28 and Equation B.29, we immediately obtain the 


expression for the inner product: 
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Qn? ntm)ln(n + 1)(4n? +4n4+3 
k(2n + 1)? (n —m!)! (2n — 1)(2n + 3) 
Finally, using the orthogonality of the L functions (Equation 3.21) : 
20 ™ 
| do f° dO sind Etec, 

0 0 

2 
_ _4Ank FY Ti? (br) + (0 + Da (EO) baerbne' (D.44) 


(2n +1)? (n—m)! 
and the integrals given by Equation B.28 and Equation B.29 the expression for the 


inner product between the L functions assume the form: 





2n*k = (n+m)! (4n? + 4n—-1) 


(2n + 1)? (n — m)! (Qn — 1)(2n +3)" (D.45) 














When we consider the inner products for m < 0, by virtue of Equations 3.10, 
3.11 and 3.12, the inner product for the (negative) —m’s will be equal to the inner 


products for the corresponding (positive) m times a non-negative factor. Thus: 





(ban [Lynd 


(Mara|Mm) = ( 


(n+m)! 


(n —m) 


! 


! 


(n +m) 


(ea cites 


) (Mun [Mn 


(D.46) 


(Ninn Nn) (eae ey 


(n+m)! 











D.10 Evaluating Inner Products when m < 0 


When m < 0 we can evaluate the inner products from the known results for m > 0. 


The m = 0 case can be evaluated using either sets of formula. We will address this 
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issue with a specific example. By virtue of Equation 3.11 we can write: 


—m)! 
(Ex|[M—mn) = / e, - Men eikr — (1) (nm m) 


<— aM ge .. (DAT 
10d (n+m)! is : : ( ) 


Now it can be easily verified that the following identities hold true: 


fe dd eimd (c’* ce e7'#) eil¢ 


20 ' 3 , ; 
if db e7*™? (c# + eW##) elle (D.48) 
0 


ae denim (ci# e“i#) el® — (D.49) 


0 


20 ‘ £ : : 
| dé eine (ci# _ e7'#) ee 
0 


ie doe”? (ef + eWi#) ef = 4 i dbe"*™? (ci# + e##) eile (D.50) 


dpe”? (ci* — e~i#) eile — ic dpe? (ci@ — ee) eae DL) 


We have 
—t jn(kr) a (cos 8) (ci — eo"*) a , (D.52) 
and 
Min: = 5 [in jac) cos 9 EACOSE) (6H 4 gis) gnind 
—t jn(kr) aces) (c'* — e-'#) aud . (D.53) 


On multiplying both sides of the above two equations by e’*** and carrying out the 
¢-integration only, the (ci# + e~i#) term in Mn +e, will have the same sign as the 
corresponding (ci# + e7i#) term in M*,-e,, but the (ci# — e7i#) term in M,,, : ez 
will have opposite sign as compared to the (ei# — e~i#) term in M*,-e,. Thus, 
overall the M,,,, - e, integral will acquire and extra negative sign as compared to 
the M*,, . e, integral after the ¢-integration is carried out. The remaining # and r 
integrations will yeild identical values in the two case, except that the e~"* factors are 
ik-r 


replaced by e“%. The second [series of e*** is used when the integral is carried out 


with e’”® instead of e~'”?, | being always greater than zero. However, the M*,,, - e, 
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integral above is precisely the one used in the evaluation of (E,|Minn). So we can 


readily obtain the inner product for negative m’s in this case by following the recipe: 


























—m)! i —i(m£1)8 
(E,|M_mn) = (-1)""1 (n m)! x ieee with | (D.54) 
(n +m)! factors changed to e(™+08, 
Similarly, since: 
i In(kr) . , ; 
Nie. = 5 n(n + 1)2 er) sin 8 P’"(cos 6) (c# + e7i#) em (D.55) 
1 d(rjn(kr) dP™ (cos 9) ia 2i\" Samed 
te Oe +e Je 
1 d(rjn(kr) Pe (cos 8) id -~id\ .-imd 
ep dr sin 9 (¢ — r : 
and 
1 (kr). sper 2 
Ninn @r = > n(n + 12 er) sin 8 P”"(cos 6) (ei# + erm) ene (D.56) 
1 d(rjn(kr) dP™ (cos 8) ib aN Gd 
ae ae Co eae ie (c +e Je 
1 d(rjn(kr) Pi (cos 8) id ~id\ ,imd 
ae ar: sin @ (e aii : ; 


now the @ integration will cause a change in the sign in the third term, but no change 
of sign for the first and second terms. In effect there will be no extra negative sign 
as that obtained in the previous case. Thus we can immediately write down the 


short-cut to obtaining the negative m inner products in this case: 





(n—m)! (Ez|Nmn) with e7*m+1)6 


E,|N—mn) = (—1)” 
Sa ae (n+m)! — [factors changed to e(™+)8, 


(D.57) 











Exactly the same considerations are valid for evaluating (E,|L_m,) since the ¢ 
parts in the integrands are identical to that used in evaluating (Ez|N—mn). So we 


have the following expression: 
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(E2|L—mn) = (—1) 


ohn yt i -i(m4£1)8 fac- 
m(z—m)! g heen with e a (D.58) 


(n+m)! — [tors changed to e(™*08, 











Now we consider the e, terms. We have: 


1 ne 6 ; 
Minne, = 5 jal) cos 8 Ey (eos?) (ci# = eW##) ee 


sin 9 


dPi"(cos@) 7 is eee 
— jr(kr) ar ie (c +e ) € , (D.59) 
and 
Me, ey = ; [=n nt cos 6 Pa(cos 8) (c# — a?) ee 
; dPi"(cos8) ¢ ig | ~id\ .-imé 
— jn(kr) a (c +e ) € . (D.60) 


So now the ¢-integration will change the sign of the first term and not change the 
sign of the second term, so that there will be no relative sign difference in the above 


two terms. Thus: 




















(n — m)! ie [Minn) with ened 
E,|M_imn) = (—1)”"+————_. x # D.61 
(E,| Poe) (n+m)! factors changed to el(m£1) 8 | ( ) 
Similarly, since we have: 
* 1 . jn(kr) . m id -i¢\ ~-imd 
IE eye ve 5 —in(n + 1) = sin 8 P™ (cos 6) (c —e ) € 
1 d(rjn(kr) dP™ (cos 9) ié ee 
at nee —e Je 
1 d(rjn(kr) Pr'(cos 9) ( i¢ -i\ ,-imd 
—im Fa Tp ae (c +e ) € ,(D.62) 


and 





1 ny: | eo 
Neva eee = 5 [inn + De en sin 6 P™ (cos 8) (c# — e7'#) ene 





1 d(rjn(kr) dP” (cos 8) id -i¢\ ime 
-iE Eon Be a Je 

. 1 d(rjn(kr) P™(cos 8) id -i¢\ imd 
+ ire Tr smd (c +e ) € (D.63) 
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So in this case, the ¢-integration will change the sign of the first two terms and not 


change the sign of the third term. The result is an overall factor of -1. Thus: 





(E,|Nimn) with e~i(m+1)6 
i(m+1)8 


factors changed to e’ ; 


(n—m)! 


(Ey|N—mn) = (ape"" (n ai m)! 


(D.64) 











Also, as discussed before, the (E,|L_,,,) integral will have the same sign inversions 


as the (E,|L_,,,,) integral and so, 














—m)! (E,|Linn) with e~*™+08 fac- 
E, |Lomn) = (1, | Balen | D.65 
FI pay (n+m)! — [tors changed to e(™08, ) 
Finally, we take a look at the e, integrals for negative m’s. We have: 
e,-M*,, = (-sin@) (—imnin( br) ECO) -om (D.66) 
sin 6 
ar P™(cos@) ; 
e (— sin 8) (in (kr) Ti (D.67) 


It is easily seen that the two expression differ in sign. So, we can immediately write 
down the transformation rule for the inner product in this case (also, since 1 = +m 


are the allowed terms, the exponential factors are e+’? ): 





mine aera an ae with e~imé al 


al (D.68) 


tors changed to e'”9. 











Both, the L and N functions have their zm factor in the ¢ components only. Sinze 
e. does not have any eg component, therefore there is not extra phase change in the 


case of L and N functions, and therefore: 





(n—m)! 


(E.|N_mn) = (ahs (n ate m)! 





(Ey|Ninn) with e~""? factors 
changed to e’”?, 








(D.69) 








(E.|L—mn) 


a 


(n—m)! 


(n+m)! 


(E,|Linn) with e~’”? factors 
changed to e’”9, 
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(D.70) 
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Appendix E 


Scattering Cross-Section and Total Cross-Section 


In this appendix we derive the formula for the scattering cross-section and the total 
cross-section of a spherically symmetric particle, subjected to plane electromagnetic 
waves. The total field outside the particle is the sum of the incident and the scattered 


fields: 


E 
H 


E, +E, (E.1) 
H; + H, (E.2) 


where the subscripts 2 and s denote incident and scattered respectively. 


We begin with the definition of the complex Poynting vector 


1 
sS 5B x H* (asterix ‘x’ indicate complex conjugate). (E.3) 


Now, since 
Eo = Ei + Eso; Eg = Eig + Eng; Ho = Hie + Heo; He = Hig + Hoo. 

So the radial component of the total complex flow vector is: 
Se 


r 


[Eo H3 — Ey Hj] 


Le el NO ee SO 


(Eve + E,a) (Ay, af: H:, | — (Eig + Ese) (Hig + H’,)| 
(HH, = FisHy) rE ; (Bo H%, 3 EH) 
es ; (Lio H*, + Eo H%, — EigH*y — EspHi) (E.4) 


Assume a spherical surface of radius R, concentric with the scatterer and enclosing 
it. The real part of S* integrated over this sphere is equal to the net flow of energy 


across its surface. Total energy absorbed by the sphere is given by 


wT 27 
W, = —Re | dé sind f ddr? S* (E.5) 
0 0 
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The first term in the right hand side of Equation E.4 measures the flow of energy 
in the incident wave. When integrated over a closed surface, this gives zero, as long 
as o2 = 0 (ie. surrounding medium is non-absorbing). The second term measures 
the outward flow of the energy scattered from the particle, and so the total scattered 


energy is 


1 Tw : 27 ‘ 
= 5 Re it dO sin 6 | dpr? (EyoH% — EspH%) (E.6) 


Since the total energy must be conserved, therefore the third term in Equation E.4 
must equal in magnitude to the sum of the absorbed and scattered energies. So the 
total energy derived from the primary wave and dissipated as heat and scattered 


radiation is given by 


Ww. = Wi+W, 
1 wT 20 
—5 Re | dO sin 0 | dd r? (EH, + EH, — BigH%g — EspHiy)(E-1) 
0 0 


Since we are considering a spherically symmetric scatterer, therefore we can choose 
the incident wave to be directed along any any arbitrary direction. We choose the 
direction of approach to the the z-axis, so that we can write the incident fields in the 


following form: 


Bi = EY {ality + iin} (8) 


n=1 


and since V xX m= kn and V x n = km, therefore, 








He’ = Vx &E; 
ao fla 
ko a : 
_ 'b§m (2)3 
on at m9 - ain ane} 





- = 


er (219-7 Fo (2) 
avi ei a E. 
Hea nol n(n 1) t! :. Si ae vin ( ") 


Similarly, the scattered fields can be written as: 





| 5 ———— 4a’ ben E.1 

s a n(n an 1) {a* mi" +2 cee ( 0) 
he (2 

H, = —>-- ae 1) fiaghon Ot tanh (E.11) 


0 [2 n=] n(n a 1) 


So we have: 


eH 1 
He 5 2”(2n + 1) {atm teh i beh 2) eo} 


a n(n +1) n oln eln 
Raa a (2 1 Pp 0 
n=1 
1d 
per : 
a ade do” 


Similarly: 


oe an(2 1 d 
Ea = ye al CineeD) fest . (—A (kr) . agen (cos 6) - sin d 


ere Ae ee 
tabs" (=f sterner) : 


sin 8 


kg (2 1 
Ae = = aCe {bs'm 7&9 — tat ; es} 
bo fig n=1 n(n oi 1) 








Wir <4 n(n +1) sin 0 


P}(cos 9) 


ky QR iM(Qn 41) {oe (—n( (ery) acts 6) ae 


162 


(rh (kr)) Spies 6) cos ‘| (E.12) 


singh, (E.13) 


iat (7 + (rhi(kr)) GPilooss) sina} (E.14) 








aE oP 
eae ph (_pO(kr)) . — P!(cos 8) - 
: wpa <= n(n + 1) ibs ( oA r)) 7 (cos #) - sin d 
; 1 d P}(cos 9) 
apa") = (ehh EN 
oe (;- dr (r nf r)) sin @ 
Therefore, 
Baty = SAD [atten ee eo 
wpe | n(n +1) con) 


L vd 


+ib8" — — (rh (kr)) - SP(cos @) - cos a1 x 


” kr dr 





s (—2)"(2n' +1) {iM sb) : SPi(cos 0)-cos¢ 


rear ni(n' +1) 


sin 6 


7 KS. 1 d 
+iats" (Fe plrentr)) 


* P1(cos 8) 


cosh (E.15) 


(E.16) 
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7 ho a ie : 4 @ 4 . 
Eg: Hig = a b> bani, fash AM(br) Goin (cos 9) : sing (E.17) 
OG een ee estes cae 
Rae ( a ha (k ») sin 9 ad lee 
(al (Ont 1 P! (cos 6 : 
ER eco 


. d * dd 
+iazs" (= ae Streuli) ag Pu (cos 6) + sin 1 : 


The ¢ integration over [0,27] gives 











20 20 

jh dé cos” d or dé sin’ d (E.18) 
0 0 
which evaluate to 7. The 6-integrations involve: 
df) APs, 1 
a d@ sin 8 (5 9 + an? 6 - P} P3) = 0 (Ee n’) (E.19) 
7 2n(n + 1) (n+1)! ie) 
2n+1 n—1)! 
7 Pe OBS = Pas ae 1 pl 

[ sino (25: r +B) = [oS (ees) (E.20) 


= Pi (cos0).P) (cos @)|" 
= 0 (since P}(4+1) =0). 


Thus we have 


| ag (EoH%y — EsgH’) = (E.21) 
whe © 2 [i"(—2)" (Qn + 1)(2n' +1) 
sae yD | n(n + 1)n'(n! + 1) e 


W fla n=1nl=1 





P}(cos 0) d 
__ sh pxsh | p(1) : (1)* SST, 1 
t( ane Or RLY (kr) «hy (kr) mane a — P2,(cos 6) 





ld : Pi(cos@) Pi, (cos 6) 
sh gxsh p,(1) ~ | 7p) n eR 
Tian dat hin '(kr) (z ap (ir) sin 8 sin 6 
1 
—ibs" pee (Geter) AGD* (kr )- opie 6): a apt (cos 9) 


kr dr i do” do” 
Ps at) 


sin 6 


ld led “od 
=pant AW (k — (rp (k — P!(cos 6 
ay, (eee (Ar))} | oa (rhe (kr) | ap Fn (cos 8) 
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P2,(cos 9) 


i ° 
i. sin 6 


= («e tc AW (kr) ali (br) © PM cos 0) - 





1 d “ad d 
—iash ats AD (kr) - (Ferien) ) apt n (cos 0) - 76 — P},(cos 8) 
ld Pi(cos@) P},(cos 8) 
psh psh rh Vik . pb be gc” nega 
TOR On (= mn ee) we (Ar) sin 8 sin 6 


sin 6 


sumac (LZ (mMary)) (EA twaldyary)) 7A + pr (cos0)) 


Due to the 6-integrals just discussed, we have the ane! and the bra aef terms 


vanish when the 6-integration is carried out. The ias*a*#" and ibs" bz" terms eee the 
same r-dependence, and on 6-integration, only the n = n’ terms survive, so one of the 
infinite summation drops out and we can write ta2"a*#" = z|as"|? and ibs*b*3" = i|b3"|?. 


Thus we have 











oe do [ d@ sin 6 (EH, = E.gH%) - (£22) 
0 (0) , : . 
on u i {ith (ie jelrteter)) ¥ 
7 1pl gpl dp! 
[sino or o v) 
-wsarcen (se geactany) f'araine (GiGr +e) 


_ aay fF 9 Sei i ash? ie (Zztrnveeny)) ) 


HIP (asthe) (Gert wry) ) 


Now the asymptotic behavior of h((kr) is given by 








1 : 

AD (kr) ~ ae (E.23) 
r 
7)". 

AY (kr) — on (B.24) 
ig YR : 

AQ (kr) a ae (E.25) 


kr 
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Therefore, 


thW (kr) - (ae sqirtt Mtr) 


in(kr) (“ = elie mA) 





(2n + 1) 
ia) es? etkr (n + 1\(- i)"e ikr __ (—i)"?2n et * 
kr (2n +1) ( i . 
ee Be wt 
=" Gas Onn 
(—ayrtyayett 


~ (kr)? (kr? ae, 


Thus the total scattered energy at the sphere at radius R is 











1 wT 2a 
WwW, = she | dé sind | dé (EuoH%, _ E.gH%s) -R 
2 0 nt) 


aN R - sh|2 sh|2 
~ One esa, (E.27) 
n=1 


= ay (2n + 1)(as”|? + [b8*|?). (E.28) 


The scattering cross section is defined as the ratio of the total scattered energy 
per second to the energy density of the incident wave. For an incident plane wave of 


unit amplitude, the mean energy flow per unit area is 


EF? 1 . 
(E x H) = 3 Vale = 5 V e2h2 ( since Kyo = 1) (E.29) 


Thus the scattering cross section, Q, is 





Q.= FE LOn + (last? + [sh (E.30) 


2 n=1 











Let us now find the sum of absorbed and scattered energies, W; : 


1 wT 20 
W, = —5Re ‘ dO sin 6 i dd (EjHt, + Eu H*id — EigH" 30 — EjgH,) FP 
0 i) 
(E.31) 
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Here we have the two combination terms: 
DFE = Eig Hy and Es Hi, — EgHjp. (E.32) 


The expressions for Eig, Eig, Hig and Hjg are similar to the ones derived explicitly 
for Eo, E.g, Hyg and H,g, with as” replaced by a‘, and bs" replaced by b¢. Also 
the spherical Hankel functions are replaced by the spherical Bessel functions. Let us 


evaluate the two terms separately. 


20 
fe : d¢ (EH, — Eig’) = (E.33) 


th, SS Sear L)(2n' +1) 
n(n + 1)n!(n' + 1) 





[2 n=1ni=1 
. Pi(cos@) d 
_ it yesh | _ 7 ()* n : 1 
( a OF <4,,(kr) eh (kr) F 7; P*(cos 6) 


* P!(cos@) P1,(cos 6) 


sin 6 sin 6 





tial, a'3" j,(kr) (Fela? or) ) 


=i ph" (ge lrintir)) RO (kr) - Pi(cos 0) - Pii(cos0) 


ld... 1d “d P1,(cos 6) 
—bi at® | ——(rjn(k —— (rh) (k 2B! (coe) 
n Gn (;- Fae ( ")) (F ra w (kr) ap in (cos ) sin 8 

1 
= («: ote 3 (kr) nl*(kr) 5 P2(cos 0) - Polos) 


: 1 * 
—ial, aXe j,.(kr) - (ee EerHsPer)) SP(cos 0) - SPA (cos 6) 


Pi(cos@) Pi,(cos 6) 





1 
SEY a bite (a ulin) -AY*(kr) - 


kr dr sin 8 sin 9 
To dene 4 ld : Pi(cos#) d 
4 «sh (1) a eet ‘ 1 
at Oe ays (; fp(rialtr))) (z Flt Pn! (br) sin @ sp Pitcos®) \ i 


Clearly, the only surviving terms after 9-integration will will be the a‘,as? and b' b3? 


terms with only n = n’ surviving. 
20 T . ' : 
= | dd i. dO sin (Ei H%, — EigH’) 
ky & . 1 ; 
TS 202m 41) {anon (iinltr) (Zr) | 


Whig 
+ bibs (inr¢an (Eatrintiry) ) | (E.34) 
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Now we consider the second term, 


20 
= | dé (EH, — ExsH) = (E.35) 


thy SS poe + 1)(2n' +1) 3 
n(n + 1)n'(n’ + 1) 





Ww fla n=1n/=1 


Pi(cos6) d 


__ sh pei, p() _, STUN Sve Dd 
i( a’ On, AM (kr) -jp(kr) smd Fat ni(cos 8) 





: bods... i Pi(cos@) P),(cos 9) 
sh _*i (1) are ; n on 
HB ei le A) (F dr (rn (ir) sin 6 sin 6 
_f{ 1 
—ibs* pri (Fe eeriter)) 7%, (kr) - SPcos 6). “Pi(cos) 


sin 8 


fled ld. og P2,(cos 9) 
_ psh xt (1) : 1 n 
Bae (= Gp hin (ir) (= ap tsn (ir) pin 08 pyro) 


. d Pi (cos 6) 
7 sh pxi p(1) * “ pil on 
(« br, AM (Ar) jn( kr) ae P, (cos 8) err ae 
1 : 
—ias* an AY (kr) : (ge riot) SPMcos 6) - TPA(cos 6) 
Pi(cos@) Pt,(cos 8) 


sin 8 sin 9 


fi 
ipso (EL tonl(br))) + Jahr) 





Poe ee) 1 0 “ Pl(cos@)  d 
+62" aX (Fe eriter)) (riot) Fa(cos 6) ap tn (oos 0) \ 


So the second integral is 
20 wT 

: dé fi dO sin @ (EoH?, — EsgH%) 
0 0 


Th > 2(2n +1) fare (ina (Zatriatin))) ) 


WH? n=1 


L, 


+ (iit (ZZtrnerery) ) | (E.36) 


Now let us consider the asymptotic behavior of the spherical Bessel functions: 





1 1 
jal kr) & = 008 (kr = 7S nr) 








r 

1 eilkr— 2H*7) 4 en i{kr- fr) 
kr 2 

1 ikr it) | ike ite 
~ ee 2 +e ‘Ee 2 


Qkr 
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f is 
~ are {eft : (—2)"*1 ae ee alee (E.37) 
Similarly, 
. ‘ thr “\n —ikry s\n 
jn-i(kr) & ee (—)” + e7*"(2) }, (E.38) 
. i thr “\n —ikry s\n 
jnti(kr) & See (—2)"+e k (2) e (E.39) 
Therefore, 
Tey, 2 1 ; 
kp ap Panlkr)) = Qn 41) {(m + 1)jn-1(kr) — rngi(kr)} 
1 “\n tkr ‘\n  —tkr 
~ aa EP + (2)"e my (E.40) 
Similarly, 
Lee) | (1) (1) 
bp dp hen (kr)) = Qn +1) {(n + 1)h,-i(kr) — nhl (kr) } 
grew thr 
~ ip (E.41) 
Therefore, 
iintkr) | EA (n ny] = sri apt bee) 
i kr dr~ ™ 2(kr)? 
1 —2ikr “\2n 
Sere ee 
1 —i2kr 
— Tere + (1) ter}. (E.42) 
Similarly, 
* (—])etl : 
th (kr) . lz . Frnt) = a 2 etkr é = , {()"e-**" + (—2)"e"*"} 
1 : : n t2kr 
= pe HH-9 -(e™ 
1 : 
= 5p [1 + (—1)"e**"]. (E.43) 
So, adding f; and fj, we obtain: 
Tk 1 ies 
——— . 2 1 
eres wpte (kr)? | a 
eee ue (—1)"+1e7#*r) +4 a eas Ol +4 (—1)"e7#*r) 
+aiath(1 + (—1)"e"*") + bitbst(1 + (—1)**4el4") ] (E44) 


t xsh 
If we let a, a**" =a 


bt bes = band 


e) nm-n 


within the square brackets as: 


Since, 


we obtain 


WwW, = 


a(l—c) + b(1+c¢) + a(1+ec*) + VI - cc) 
(a+a*)+(b4+8)+cb-—a)—c(h —a*) 
2Re(a) + 2Re(b) + 22Im(c(b— a)). 


R? 
W, = > Relh + 1] 5 


a €2 = i ks i L*S 
ara rae 200 +1) (a‘,a% A bt Be *) : 


and since the total cross section is defined as 


therefore: 











Q: 


_ 2x 


= t ox«sh a pxsh 
2 Re S/(2n +1) (a‘,a% + bb ) 


n=1 
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(—1)"e7*" = ¢, then we can write the expression 


(E.49) 


(E.50) 
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Appendix F 


A Vector Identity 


We give the proof that the following vector identity holds good in the spherical polar 


coordinate system: 





VxVxA=V(V-A)-V7A (F.1) 











In spherical polar coordinates, (r, @, ¢), unit vectors at point P are e,, eg and 


e,. Unlike the unit vectors in rectangular coordinates (e,, e,, €.), €-, €9, €4 are not 


constant vectors: 


e, 
e, 


e, 


Similarly, 


Thus, 


(e, - @,)e, + (e, -e,)e, + (e, - ez )e, 


sin@ cos¢ e, + sind sing e, + cosé e, 








( sin® cos, sind sing, cosé ) (Cartesian components) (F.2) 
es = (cosé cos¢, cosé sing, —siné ) and (F.3) 
es = ( —sind,cos¢,0) (F.4) 

de, es Jeg 

Or Or Or 

a = ( cos@ cos¢, cos@ sing, —sin# ) = eg 
de, gt at 

— =( —sin@sin¢, sin® cos¢, 0 ) = sin? eg 
a6 

es 


=( —sin@ cos¢, —sin@ sing, —cos# ) = —e, 


00 
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= = ( —cos@ sind, cos@ cos, 0 ) = cosé eg 
Jeg 

—_ = =O 

0e ( 0, 0, 0 ) 

Jeg ‘ : 

ob = ( —cos¢, —sing, 0 ) = —sin# e, — cosé eg 


Let us define an arbitrary vector function A: 
A = Adt, g, d)je, + Aa(r, 0, d)es + Ag(r, 0, deg. 


The Laplacian operator V? in spherical polar coordinates is given by: 


ee eer mee feces ne 
= Pisind |ar \ ar} | OO \” a6) | Od \sind Od 


The first term(operator) does not act on e, , eg or eg by virtue of Eqn F.5. If we carry 


2 











out straight forward differentiation, we obtain: 








< (resine < (A,e, } = \¥ (r2sne =) he. (F.5) 
< (resin © (Aves ) 2 iz (r2sina oY be (F.6) 
< (r2sine < (Ages \ — iz (r2sne “se es (F.7) 
Similarly for the @ and ¢ parts, we obtain: 
< (sino < (A,e, ) = -sind A,e, + sind oe 
ie ee (sind a} ey + ee (sine a e, (F.8) 
4 (sno 4 Be ) = a (sind Ay) e, — sind Ages 


2s 2 OAg Of. ,0A¢ 
—sin§ oT + a (sine ) eg (F.9) 


< (sine < (Ages \ = 135 (sine oe es (F.10) 


Li2 

















- (= = (Are, ) = -sin A,e, —cosd A,ey 
+ (=a i (Aves ) = —cosé Ase, — cot Aces 
1 tes +155 (aera) ber Fa) 
i" (saa 3g (Ae ) se Oe <0 Aves - 2 Fe, 








cos? JAg 0 1 OAg 
ay ae ss in (= Oe )hee ee) 


Thus we observe that V7A = V?(A,e, + Ageg + Age ) will have components 
along e, , egand eg. Therefore (V?A)-e, will be given by the sum of all the 
e, components obtained above (Eqns F.5,F.8,F.11): 


OAs ~20A; cos6 OA, 10°A, 2 
































VA) ea = — —A, 
( )se Or? 2 r Or - r2sin0 O60 - r? 06? r2 
2 2 1 ae 2 
a OAg : goed ae O7A 7 OAg (F.14) 
r? 06 r*sin@ r?sin*@ 06? ~—r*sin@ O@ 
Similarly, 
2 OA, 2 OAg 0? Ag Ag cos6 OAg 
V7A)> = = 
( yeee ¢? ~O0 zn r Or - Or? r2—s r*gsin@ = O60 
1 0? Ag cos?6 1 OAs 2cos OAs 
— — wl 
ue 002 r2sin26° * a r’sin?@ 062 —r?sin?@ 06’ ey 
and 
2 OA, 2cos6 OAg 2 0Ag 0? Ag 
V°A)- = 
( )-e¢ rsind 06 ~ r*sin?6@ 06 — r Or Or? 
cos@ OAg 1 OP Ag 1 OP Ag Ag 
— PL 
r2sin?@ 06 r? 06? 7 r2sin?6 Od? r2sin?6 (P16) 
Now, since 
OA, 2 cos@ 1 OAg 1 OAs 
-A= A, : ; ; £. 
” Or x r 3 rsind © z r OO  rsin@ 0¢ ue) 
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Therefore, [V(V- A)]-e, = 2(V- A) is given by: 









































OPA, 2 20A, cos6 cosO OAg 
. A . r — mag r —— " ° 
Nene ie Or? vag r Or rsind © | rsind Or 
1 0Ag 10°Ag 1 OAs 1 @PAg (F.18) 
r2 00 rdOr0@ r*sin@ O¢6 — rsin? 0¢0r 
Similarly, [V(V-A)]-e9 = +5(V- A) is given by: 
1 0° A, 2 OA, cos6 OAg 1 Ag 
V(V-A)]- = = 
IV( ee r000r  r? 00 r*sin@ 060 — sin?@ r? 
1@A 2A 06 OA 
OP Ag ! Ag _ cos O 6 (F.19) 
r? 06? r?sin6 0006 r?sin*@ 0¢ 
and, [V(V-A)]-e, = vane ag (V - A) is given by: 
1 O?A, 220A, cosO OAg 1 07 As 
Sa ES rsin0 So r Oo rsin8@ 06 rd¢08 
O° Ag 
— F.2 
rsind Od? | 20) 
Now, 
7 10Ag | cosé 1 OAs 
Mee |? 00 * Fane ° rsind mle. 
1 OA, Ag OAS 
z 5 (ero) op Or | mt 
OAg Ag 104A, 
- = ror vaca et) 
Therefore, 
107Ag cos OAg 1 OAs cosé lo? A: 
A]-e, : : = 
EI E O00r  rsin@ Or r2 00 | rsind © 2? 06? 
2 2 
2 cone OA, z : O° A, i : OAg z OP Ag (F.22) 
r’sind 00 ~~ r?sin?@ 0¢2_ ——sr?sin6 Od — rsinO Od0r 
and, 
1 OP Ay cos@ OAg 1 0? Ag 1 OAg 
Al. = Z 
Pex Vie ces sana 0600 : r2sin’?@ O06 ~~ r2sin?@ 062s: Or 
O?Apg 10Ag 107A, 
— — F.2 
Or? r Or r Or0o 2s) 
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and 


1 @A, 204A, GA, 1 O?Ag 
rsind Ordd or Or Ar??? 





[Vx VxA]-e, = | 


cos6 OAg Ag cos OAg 1 6?As 
=a aa +— (F.24) 
r2sind OO r*sin°@ ~~ r2sin*O 06 — r?sin@ 000d 





Therefore we find, by comparing term by term, that: 


[Vx VxA]-e, = [V(V-A)-V?A]-e, 
[VxVxA]-eg = [V(V-A)-V?A] -e5 
[Vx VxA]-eg = [V(V-A)-V?A] -e, 


Hence the proof. 


